Graph Ricci curvatures reveal atypical functional connectivity in autism spectrum disorder

While standard graph-theoretic measures have been widely used to characterize atypical resting-state functional connectivity in autism spectrum disorder (ASD), geometry-inspired network measures have not been applied. In this study, we apply Forman–Ricci and Ollivier–Ricci curvatures to compare networks of ASD and typically developing individuals (N = 1112) from the Autism Brain Imaging Data Exchange I (ABIDE-I) dataset. We find brain-wide and region-specific ASD-related differences for both Forman–Ricci and Ollivier–Ricci curvatures, with region-specific differences concentrated in Default Mode, Somatomotor and Ventral Attention networks for Forman–Ricci curvature. We use meta-analysis decoding to demonstrate that brain regions with curvature differences are associated to those cognitive domains known to be impaired in ASD. Further, we show that brain regions with curvature differences overlap with those brain regions whose non-invasive stimulation improves ASD-related symptoms. These results suggest the utility of graph Ricci curvatures in characterizing atypical connectivity of clinically relevant regions in ASD and other neurodevelopmental disorders.

www.nature.com/scientificreports/ Graph theory and network analysis provide objective, data-driven measures to analyze the topological architecture and connectivity patterns (human 'connectome') in the human brain 14,[19][20][21][22][23] , and can provide us with deeper insights about the functional, structural and causal organization of the brain 23 . Remarkably, many previous studies on ASD have utilized graph-theoretic analysis of rs-fMRI functional connectivity networks (FCNs) to differentiate ASD from typical development 18,[24][25][26][27][28][29][30][31][32] , and furthermore, some of these studies 18,24,30,31 made use of the ABIDE-I dataset. These studies investigated network characteristics such as small-worldness, modularity, clustering, efficiency, rich club organization and connection densities of the FCNs in ASD versus typical development, and reported atypical functional organization in ASD both globally and at the level of individual brain regions.
In recent years there has been an increasing interest in the development of geometric tools for analyzing complex networks 33 , which enables the study of higher-order correlations in networks beyond pairwise interactions [34][35][36] . A fundamental concept in geometry is Ricci curvature 37 , which quantifies the extent to which a space differs from being flat. Various nonequivalent definitions of graph Ricci curvature have been proposed [38][39][40][41][42] with an aim to capture the key properties of the classical Ricci curvature. Different notions of graph Ricci curvature have found applications in diverse areas, such as differentiating gene co-expression networks of cancer cells and healthy cells 43 , identifying crashes and bubbles in financial networks 44,45 , and community detection in complex networks 46,47 . Ollivier-Ricci curvature (ORC) 40 and Forman-Ricci curvature (FRC) 39,41 are two widelyused notions of graph Ricci curvature.
Notably, graph Ricci curvatures have also been applied to structural and functional connectivity networks of the human brain. Farooq et al. 48 applied ORC to brain structural connectivity networks to identify robust and fragile brain regions in healthy subjects. They also show that ORC can be used to identify changes in brain structural connectivity related to ASD and healthy aging. Simhal et al. 49 used ORC to measure changes in brain structural connectivity of individuals with ASD before and after the infusion of autologous umbilical cord blood. ORC has also been used to study differences in brain structural connectivity networks of cognitively impaired and non-impaired multiple sclerosis patients 50 . Recently, Chatterjee et al. 51 used a version of FRC to determine the changes in brain functional connectivity related to attention deficit hyperactivity disorder (ADHD). Additionally, FRC has been used to analyze task-based fMRI data 52 as well as to predict the intelligence of healthy human subjects 53 . Most of these studies have also contrasted graph Ricci curvatures with standard network measures such as clustering coefficient and node betweenness centrality, and showed that graph Ricci curvatures can provide new information about brain connectivity organization. However, a systematic evaluation of the ability of graph Ricci curvatures to characterize atypical brain functional connectivity in ASD and other neurodevelopmental disorders is lacking.
In the present work, we expand the scope of curvature-based analysis for characterizing brain connectivity, by systematically applying graph Ricci curvatures to study atypical functional connectivity network organization in ASD. For this purpose, we utilized raw resting-state fMRI images of 1112 subjects from the ABIDE-I dataset and obtained FCNs for each subject by implementing a uniform preprocessing pipeline and thorough quality assessment (QA) checks. We employ FRC and ORC to compare the FCNs of individuals with ASD relative to typically developing individuals (TD), and evaluate the role of these curvature measures as indicators of atypical functional connectivity in ASD. We analyzed the brain-wide changes in FCNs by comparing average edge curvatures across the two groups, and analyzed the region-specific changes in FCNs by comparing node curvatures across the two groups. Then, we used meta-analysis decoding with respect to a large database of fMRI studies, to determine if those regions showing curvature differences are also associated to those cognitive domains known to be impaired in ASD, e.g. social cognition. Finally, we determined if those regions showing curvature differences overlapped with those regions whose non-invasive stimulation with transcranial magnetic stimulation (TMS) 54 and transcranial direct current stimulation (tDCS) 55 are reported in the literature to result in improvement of ASD-related symptoms.

Results
The primary goal of this study is to evaluate the utility of two notions of graph Ricci curvature, namely Forman-Ricci curvature (FRC) and Ollivier-Ricci curvature (ORC), that have been recently ported to the domain of complex networks, as indicators of atypical topological organization in resting state FCNs of individuals with ASD. For this purpose, we analyzed spatially and temporally preprocessed rs-fMRI images of 395 individuals with ASD and 425 TD individuals from the ABIDE-I dataset as described in "Methods". The demographic and clinical information for these subjects is summarized in Table 1. Figure 1 is a schematic summarizing the processing pipeline for rs-fMRI data used in this study. Furthermore, Supplementary Table S1 gives detailed information on the quality assessment and exclusion criteria for rs-fMRI dataset. Subsequently, 200 regions of interest (ROIs) or nodes were defined in the brain using the Schaefer atlas, and a 200 × 200 functional connectivity (FC) matrix was generated for each subject by computing the Pearson correlation coefficient between the time-series of all pairs of nodes. Thereafter, by combining maximum spanning tree (MST) and sparsity-based thresholding, we constructed FCNs over a wide range of graph densities between 0.02 or 2% edges and 0.5 or 50% edges, with an increment of 0.01 or 1% edges (see "Methods"). In a nutshell, we generated and analyzed 49 FCNs for each of the 820 subjects in the ABIDE-I dataset considered in this study.
Brain-wide changes in functional connectivity networks. To investigate the differences in the global organization of FCNs between the ASD and TD groups, we computed the average edge FRC and average edge ORC across the 49 FCNs across the graph densities 2-50% for each subject. To compare the average edge curvatures at each graph density between the ASD and TD groups, we employed a two-tailed two-sample t test followed by FDR correction (see "Methods"). In Fig. 2a,b, we show the differences in average edge FRC and average edge ORC, respectively, between the ASD and TD groups across the graph densities 2-50%. We find that www.nature.com/scientificreports/ average edge FRC is significantly lower ( p < 0.05 , FDR-corrected) in the ASD group compared to the TD group in the graph density range 5-50% (Fig. 2a). Similarly, we find that average edge ORC is lower ( p < 0.05 , FDRcorrected) in the ASD group compared to the TD group albeit the differences were insignificant ( p > 0.05 , FDRcorrected) in the graph density ranges 2-5% and 24-33% (Fig. 2b). Although the directionality of the differences with the two discrete Ricci curvatures is the same for the two groups, that is, average edge curvature in the ASD group is lower than that in the TD group, it is important to emphasize that the two discrete Ricci curvatures capture different aspects of the classical Ricci curvature, and thus, cannot serve as alternative measures across different types of networks. Specifically, ORC captures the volume growth property of the classical Ricci curvature whereas FRC captures the geodesic dispersal property 42 . While ORC has a deeper correspondence with the classical Ricci curvature, FRC is based on a simple combinatorial expression which is significantly faster to compute in larger networks. After comparing the average edge curvatures between FCNs of ASD and TD groups, we find that the statistical test (t test followed by FDR correction) yielded lower p-values after FDR correction for average edge FRC compared to average edge ORC across most of the considered graph densities (Supplementary  Table S2). In other words, the differences between FCNs for the two groups are more pronounced for the average edge FRC than average edge ORC.
To gain a deeper understanding of the altered global organization of FCNs between the ASD and TD groups, we also compared six other global network measures, namely, average clustering coefficient, modularity, average shortest path length, average node betweenness centrality, global efficiency and average local efficiency. We find that the average clustering coefficient is significantly lower ( p < 0.05 , FDR-corrected) in the ASD group compared to the TD group in the graph density range 2-50% (Fig. 2c). Moreover, our results for clustering coefficient are consistent with results from previous studies that have employed graph-theoretic measures to analyze resting state FCNs in ASD 26,28,31 . Thereafter, we find that the modularity of the FCNs is significantly reduced in the ASD group compared to the TD group in the graph density range 2-50% (Fig. 2d), and our results are consistent with the results from previous studies 26,31 . Further, we find that the average shortest path length of the FCNs is significantly lower in the ASD group compared to the TD group in the graph density range 5-31% (Fig. 2e), and our results are consistent with results from previous studies 26,28 . Lastly, we find that average node betweenness centrality is significantly lower in the ASD group compared to the TD group in the density range 5-31% (Fig. 2f).
Furthermore, we have computed two global measures that characterize how efficiently information is exchanged within a network, namely global efficiency and average local efficiency. We find that global efficiency is significantly higher (p < 0.05, FDR-corrected) in the ASD group compared to the TD group in the graph density range 4-31% ( Supplementary Fig. S1). Note that the direction of the effects observed for global efficiency is opposite to the direction of effects observed for average shortest path length (Fig. 2e), as global efficiency is defined as the average of reciprocal shortest path lengths between all pairs of nodes in a network. Moreover, our results for global efficiency are consistent with the results from previous studies 26,28,31 . We find that average local efficiency is significantly lower in the ASD group compared to the TD group in the graph density range 2-50% ( Supplementary Fig. S1). Note that the results for average local efficiency are similar to the results of average clustering coefficient (Fig. 2c), since the two network measures are closely related to each other. Moreover, our results for average local efficiency are consistent with results from previous studies 26,28,31 . Table 1. Summary of demographic and clinical information for the 820 subjects from ABIDE-I project that fulfil the inclusion criteria and were selected for network analysis in this study. 395 subjects belong to the autism spectrum disorder (ASD) group and 425 subjects belong to the typically developing (TD) group. The subjects in both groups are age-matched ( p = 0.835 ). Handedness data were missing for 137 ASD and 148 TD subjects. ADI-R social data were missing for 120 ASD participants. ADI-R verbal data were missing for 119 ASD participants.  Figure 1. Schematic diagram summarizing the rs-fMRI processing pipeline employed in this study. The raw fMRI scans undergo four steps in spatial preprocessing, namely, motion correction, slice-timing correction, outlier detection, and direct segmentation and normalization. The raw structural MRI scans are normalized to the Montreal Neurological Institute (MNI) space, and segmented into grey matter, white matter and cerebrospinal fluid (CSF) areas. In the temporal preprocessing or denoising step, the BOLD time series of each voxel is extracted and the remaining physiological and motion confounds are removed using linear regression. The confounds include white matter and CSF masks, subject-motion parameters and outlier scans. The residual BOLD time series of each voxel undergoes a high-pass filtering at 0.008 Hz. The Schaefer atlas is used to parcellate the brain into 200 regions of interest (ROIs) and the mean time series for each ROI is computed. Finally, Pearson correlation coefficient is computed between all pairs of ROIs, resulting in a 200 × 200 functional connectivity (FC) matrix.
Thorough quality assessment (QA) checks were implemented both before and after preprocessing. In this figure, the head icon under denoising section is made by Freepik from flaticon.com (https:// www. flati con. com/ autho rs/ freep ik). www.nature.com/scientificreports/ Region-specific changes in functional connectivity networks. Given the significant differences in FRC and ORC of the entire brain between the ASD group and the TD group, we evaluated node-level curvature differences in the FCNs, and determined how these differences are distributed across the 7 resting state networks (RSNs) in the brain. For this purpose, we first computed node FRC and node ORC for all the 200 nodes, across the 49 FCNs with graph densities 2-50% for each subject. Second, to identify the set of nodes that show significant differences between the ASD and TD groups, we compared the area under the curve (AUC) of the node FRC and  www.nature.com/scientificreports/ the node ORC for each node using a two-tailed two-sample t test followed by FDR correction (see "Methods"). In Fig. 3 and Supplementary Fig. S2, we show the nodes or regions that exhibit significant differences ( p < 0.05 , FDR-corrected) in FRC and ORC, respectively, between the ASD and TD groups. We identify 83 regions that show significant between-group differences in FRC and 14 regions that show significant between-group differences in ORC. For FRC, the significant regions are spread across the 7 RSNs. However, they are mainly concentrated within 3 RSNs namely, default network (26 significant regions), somatomotor network (30 significant regions) and salient ventral attention network (13 significant regions). In the default network, RH_Default_pCunPCC_2 Detailed information about the abbreviations of region names can be found at: https:// github. com/ Thoma sYeoL ab/ CBIG/ tree/ master/ stable_ proje cts/ brain_ parce llati on/ Schae fer20 18_ Local Global 56 . Thus, our node-level results for FRC and ORC suggest that the nodes or brain regions showing significant differences were not distributed evenly across the 7 RSNs, but concentrated within the default network, somatomotor network and salient ventral attention network. After applying both FRC and ORC to resting state FCNs in ASD and TD groups, we find that the betweengroup differences in FRC of the FCNs are more pronounced compared to ORC, both at the global-level and the node-level. Therefore, we mainly focus on the nodes identified using FRC in further analyses. We reiterate that the two discrete Ricci curvatures capture different aspects of the classical Ricci curvature, and thus, neither of the two measures can be treated as an alternative to the other.

Scientific
We also evaluated the node-level differences in two standard network measures namely, clustering coefficient and node betweenness centrality. Notably, ORC is related to clustering in networks 42 . We identify 78 brain regions that show significant differences ( p < 0.05 , FDR-corrected) in clustering coefficient, and 4 brain regions that show significant differences ( p < 0.05 , FDR-corrected) in node betweenness centrality (Supplementary Table S3). The brain regions identified by clustering coefficient are concentrated in three RSNs namely, default network, somatomotor network and salient ventral attention network ( Supplementary Fig. S3). Further, we computed the overlap between sets of significant brain regions identified by each of the four node-level network measures used in our study. First, we found 8 brain regions that are commonly identified by both FRC and ORC, . Second, we found 71 brain regions that are commonly identified by FRC and clustering coefficient. Third, we found 5 brain regions that are commonly identified by ORC and clustering coefficient. Fourth, we found that 1 brain region is commonly identified by FRC and node betweenness centrality. Fifth, we found 2 brain regions that are commonly identified by ORC and node betweenness centrality. We do not find any brain regions that are commonly identified by clustering coefficient and node betweenness centrality.
Behavioral relevance of region-specific changes using meta-analysis decoding. As discussed previously, we identify 83 brain regions that show significant between-group differences in FRC. Subsequently, we partitioned the set of significant brain regions according to their respective RSNs and determined the cognitive domains associated to the significant regions in each RSN using Neurosynth meta-analysis (see "Methods"). The Neurosynth analysis enables identifying the cognitive domains associated to the significant regions in an RSN more rigorously than just assuming it as the putative functional role of that RSN. We limited the Neurosynth analysis only to default network, somatomotor network and salient ventral attention network, since a considerable number of regions are detected in these RSNs and since these regions are nearly bilaterally symmetrical. Figure 4 shows the significant brain regions separately for each of the three RSNs and the associated word clouds highlighting the behavioral relevance of the significant regions in each RSN. Supplementary Table S4 lists the significant brain regions and the terms associated with all seven RSNs. The word cloud for default network shows terms associated with social cognition and memory (Fig. 4a). For the somatomotor network, we find terms associated with movement (Fig. 4b). For salient ventral attention network, we find terms associated with movement and language (Fig. 4c). Notably, previous studies on ASD population have found impairments in social cognition 3,57,58 , memory 3,4,59,60 , movement 5,61-63 and language [64][65][66] . Hence, we find that those regions exhibiting ASD-related curvature differences are also associated to those cognitive domains known to be impaired in ASD.
Subsequently, we determined if there is a relationship between the FRC of brain regions that showed significant differences and clinical scores for symptom severity in ASD. To do this, we related the FRC of just the brain regions which showed significant differences in each RSN with the behavioral function associated with that RSN, as determined by the Neurosynth meta-analysis decoding (see "Methods"). First, we used ADI-R social score as a measure of social cognition and related this score with the FRC of regions in the default network. Second, we used ADI-R verbal score as a measure of language and related this score with the FRC of regions in the salient ventral attention network. We did not find any nodes that showed significant correlations between FRC and clinical scores after FDR correction. Prior to FDR correction, FRC for the node LH_Default_Temp_1 (− 47, 8, − 32) in the default network was positively correlated with ADI-R social score ( r = 0.122 , p = 0.044 ), and FRC for the node LH_SalVentAttn_ParOper_1 (− 56, − 40, 20) in the salient ventral attention network was positively correlated with ADI-R verbal score. We also repeated the analysis for the brain regions with significantly different clustering coefficient values in the 2 RSNs namely, default network and salient ventral attention network. Similar to FRC, these brain regions show behavioral relevance and are associated to social cognition and memory in default network, and movement and language in salient ventral attention network ( Supplementary  Fig. S4). Thus, we correlated the clustering coefficient values of significant brain regions in default network with the ADI-R social score and the clustering coefficient values of significant brain regions in salient ventral attention network with the ADI-R verbal score. In this analysis for clustering coefficient, we did not find any significant correlations with both ADI-R social scores and ADI-R verbal scores after FDR correction. To sum up, neither FRC nor clustering coefficient show evidence for a relationship with symptom severity in individuals with ASD.

Agreement of results from node-level network analysis to TMS/tDCS literature. In addition
to the meta-analysis decoding, we performed one more analysis to determine the agreement of our results with relevant previous neuroimaging literature. Specifically, we determined the overlap between those brain regions showing FRC differences and those whose non-invasive stimulation using TMS or tDCS, resulted in improvement of ASD-related symptoms. To do this, we performed a literature search on PubMed to identify the set of brain regions whose non-invasive stimulation using TMS or tDCS yielded positive effects on ASD symptoms. The exact details of the PubMed search query are provided in Table 2. Then, we compared this set of brain regions to those with altered FRC values in resting-state fMRI FCNs of individuals with ASD. Figure 5 summarizes the workflow we employed to collect and classify the eligible articles from the literature survey.
The studies employing TMS have reported positive effects in ASD-related symptoms after stimulating 4 target regions, namely, premotor cortex, dorsolateral prefrontal cortex (DLPFC), pars triangularis, and pars opercularis. The studies employing tDCS have reported positive effects in ASD symptoms after stimulating 2 target regions, namely, DLPFC and left primary motor cortex. Supplementary Tables S5 and S6 provide a detailed summary of the TMS and tDCS studies, respectively. Note that the target regions in these experiments are cortical regions that are defined differently from the ROIs (or nodes) defined in our study that are a part of the Schaefer 200 parcels atlas 56 . Therefore, in order to compare the results of our node-level analysis with the effects of stimulating the target regions, we mapped the Brodmann areas that correspond to target regions 67,68 to the 200 Schaefer ROIs (see "Methods" and Supplementary Table S7).
Based on the data collected from previous NIBS experiments (Supplementary Tables S5 and S6), we identified five target regions that show evidence for improvement in behavioral or cognitive symptoms associated with ASD following TMS or tDCS, namely, premotor cortex, pars triangularis, pars opercularis, DLPFC and left primary motor cortex. These five target regions correspond to Brodmann areas 6, 45, 44, 9,   ROIs also show significant ASD-related differences in FRC and 13 ROIs show significant ASD-related differences in clustering coefficient. None of these 31 ROIs show significant ASD-related differences in ORC or node betweenness centrality. A visual representation of these ROIs is provided in Fig. 6. Notably, the 18 ROIs with significant ASD-related differences in FRC are a superset of the 13 ROIs with significant between-group differences in clustering coefficient. These results serve as a form of validation, based on the literature on outcomes from TMS and tDCS experiments, that FRC identifies clinically relevant brain regions underlying ASD. Further, FRC identifies some regions that might be clinically relevant in ASD, but are not identified by other node-based network measures, e.g. clustering coefficient. Table 3 lists the target regions that show improvement in clinical symptoms associated with ASD following TMS or tDCS, the corresponding ROIs in the Schaefer atlas that show significant between-group differences in node-level network measures, the network measure that captured the differences, and the experimental studies that report the effects.

Discussion
Graph Ricci curvatures have not been previously applied to study atypical resting-state functional connectivity in ASD. In the present work, we used two notions of graph Ricci curvature, namely Forman-Ricci curvature (FRC) and Ollivier-Ricci curvature (ORC) to compare the resting-state FCNs of individuals with ASD relative to TD individuals. We found that average edge curvature can effectively distinguish the whole-brain functional connectivity of individuals in the ASD and TD groups. Additionally, we studied the differences in node curvature between the two groups and identified specific regions in the brain with atypical functional connectivity in ASD. Notably, we found that brain regions with altered FRC in functional connectivity networks of individuals with ASD, were also associated in the fMRI literature to those cognitive domains known to be affected in ASD. Further, we observed an overlap between the set of regions with altered FRC in functional connectivity networks of individuals with ASD, and those brain regions whose non-invasive stimulation in TMS/tDCS experiments resulted in improvement of ASD-related symptoms. We acquired rs-fMRI scans of 1112 participants as provided by the ABIDE-I project 18 . The large sample size of the ABIDE-I dataset offers substantial statistical power, thereby increasing the reliability of the reported results 2,14,18 . We preprocessed each scan using the CONN functional connectivity toolbox 69 , implementing thorough quality assessment (QA) checks both before and after preprocessing. For each participant, we generated a 200 × 200 functional connectivity (FC) matrix using Schaefer atlas 56 and constructed FCNs with a broad range of edge densities using a maximum spanning tree (MST) followed by sparsity-based thresholding. MST-based network construction is particularly useful for network analyses since it ensures the resulting network is always connected. Similar network construction approaches involving spanning trees have previously been used for financial networks 44,45 and brain FCNs 70 .
After comparing the average edge curvatures of the FCNs in the ASD and TD groups, we found reduced average FRC and average ORC in individuals with ASD. Similar analysis using standard network measures revealed reduced average clustering coefficient, reduced modularity, reduced average path length, reduced average node betweenness centrality, increased global efficiency and reduced average local efficiency. All the standard network measures except node betweenness centrality have previously been used to study brain-wide changes in functional connectivity in ASD 26,28,31 and our results are in agreement with previous findings. However, the changes in graph Ricci curvatures have not previously been studied for FCNs in ASD. Our results illustrate the sensitivity of graph Ricci curvatures, especially FRC, in discriminating the resting state FCNs of individuals with ASD compared to TD.
After comparing the node curvatures of the FCNs in the ASD and TD groups, we identified 83 brain regions that are significantly different in FRC and 14 brain regions that are significantly different in ORC between the two groups. FRC and ORC identify 5 common regions. Moreover, we found that these regions are bilaterally symmetrical and mainly concentrated in 3 RSNs namely, default network, somatomotor network and salient ventral attention network. Previously, Farooq et al. 48 have used ORC to compare structural connectivity networks of individuals with ASD relative to TD, and showed that regions with significant difference in ORC are present in visual, dorsal attention, ventral attention areas and temporal lobe. Our results from comparing ORC of resting Table 2. Literature search for non-invasive brain stimulation studies in ASD. Detailed summary of the electronic search query on PubMed that we used to obtain our original corpus of articles that perform noninvasive brain stimulation on individuals with ASD.

Search date Search query Search filters Source
October 18, 2021

((transcranial magnetic stimulation) AND (autism)) OR ((transcranial magnetic stimulation) AND (Asperger)) OR ((transcranial magnetic stimulation) AND (PDD NOS)) OR ((transcranial direct current stimulation) AND (autism)) OR ((transcranial direct current stimulation) AND (Asperger)) OR ((transcranial direct current stimulation) AND (PDD NOS)) OR ((transcranial alternating current stimulation) AND (autism)) OR ((transcranial alternating current stimulation) AND (Asperger)) OR ((transcranial alternating current stimulation) AND (PDD NOS)) OR ((TMS) AND (autism)) OR ((TMS) AND (Asperger)) OR ((TDCS) AND (autism)) OR ((TDCS) AND (Asperger)) OR ((TACS) AND (autism)) OR ((TACS) AND (Asperger))
year : no filter, article attribute : no filter, language : no filter We performed meta-analysis decoding, based on the fMRI literature, to determine if those brain regions with altered FRC in individuals with ASD, were also associated to those cognitive domains that are known to be affected in ASD. We found that brain regions with altered FRC were associated to the cognitive domains of social cognition, memory, movement and language. Farooq et al. 48 used ORC to compare structural connectivity networks of individuals with ASD relative to TD, and showed that regions with significant difference in ORC are related to semantic memory, socially relevant memories, emotions and visual perception. Since these results were obtained with structural connectivity networks, it is difficult to compare our results to these. However, each of the cognitive domains suggested by the meta-analysis decoding are known to be affected in ASD 3-5,57-66 .  Figure 5. Summary of the workflow employed to compile data from non-invasive brain stimulation (NIBS) experiments. The workflow is presented according to PRISMA statement 99 . First, we identified 235 potential records from PubMed. Second, we filtered the articles based on title and abstract. Third, we scanned review articles for more records and looked for existing databases for additional data. After performing the above steps, we were left with 84 potential NIBS studies. Finally, we classified the studies based on the stimulation technique (TMS/tDCS/tACS) and screened the studies individually for eligibility. We were left with 19 TMS studies and 12 tDCS studies, which were used to extract experimental data. www.nature.com/scientificreports/ In addition to meta-analysis decoding based on the fMRI literature, we assessed the agreement of our results with those of TMS/tDCS experiments involving individuals with ASD. We found that brain regions with altered FRC in individuals with ASD overlap with those brain regions whose non-invasive stimulation with TMS/ tDCS have been reported to result in improvement of ASD-related symptoms. We further note that the set of regions with altered values of other node-based network measures, e.g. clustering coefficient, which overlap with those regions identified by non-invasive stimulation, are a subset of the set of regions identified by FRC. To our knowledge, this is the first instance of using results from TMS/tDCS experiments as a form of validation of results from a graph-theoretic analysis of brain functional connectivity networks. Hence, these results suggest that FRC captures atypical connectivity of clinically relevant brain regions underlying ASD. Further, FRC might capture atypical connectivity not captured by other node-level network measures such as clustering coefficient. These results commend the use of graph Ricci curvatures as a source of hypotheses about clinically relevant brain regions underlying ASD, which can then be tested by stimulating these regions with non-invasive technologies, e.g. TMS [71][72][73] .
There was a significant difference in the IQ scores between the ASD and TD groups, which introduces a potential confound to our analysis. Previous studies involving graph-theoretic analysis of rs-fMRI scans in the ABIDE dataset have matched the groups on IQ scores 30,31,74 . However, we choose to include all subjects in our analysis rather than sub-selecting ASD subjects according to IQ, since sub-selecting would make ASD cohort less representative and the results of our analyses would be less generalizable to the typical ASD population 75 . Confining the analysis to IQ-matched ASD subjects would include only high-functioning ASD subjects in the analyses, hence, our results would not be generalizable to ASD subjects whose cognitive functioning is more severely affected. In addition, we have not included IQ as a covariate while comparing the global and local network    75 have shown that using IQ as a matching variable or covariate during studies of neurodevelopmental disorders could lead to anomalous findings about neurocognitive function.
To sum up, we find that geometric notions of graph Ricci curvature can be effectively used to determine global and node-level changes in functional connectivity networks of individuals with ASD. Importantly, we present a validation based on TMS/tDCS literature, to suggest that graph Ricci curvatures, particularly FRC, are sensitive to atypical functional connectivity of clinically relevant brain regions underlying ASD. The methods used in the present work could further be applied to study functional connectivity networks in other atypical populations. Additionally, since graph Ricci curvatures are fundamentally defined on edges, future studies could be aimed at devising edge-based methods to analyze brain functional or structural connectivity.

Methods
In this section, we describe the methodology used to construct the resting-state functional connectivity networks (FCNs) of individuals with autism spectrum disorder (ASD) and typically developing (TD) individuals, from raw resting-state functional MRI (rs-fMRI) images acquired from the Autism Brain Imaging Data Exchange I (ABIDE-I) project 18 . Note that TD individuals are the healthy subjects. First, raw rs-fMRI data were spatially and temporally preprocessed using the CONN toolbox 69 . Second, we parcellated the brain into 200 regions of interest (ROIs) or nodes using the Schaefer atlas 56 and a 200 × 200 functional connectivity (FC) matrix was generated for each subject. Third, we filtered the FC matrix using a maximum spanning tree (MST) based approach followed by sparsity-based thresholding to construct FCNs for each subject.
Participants and imaging dataset. From the ABIDE-I project 18 , we obtained raw rs-fMRI and anatomical data for 1112 participants (age range = 7-64 years, median = 14.7 years), comprising 539 individuals with ASD and 573 age-matched TD individuals. ABIDE-I project is an international effort by 17 imaging sites that have collectively shared rs-fMRI, anatomical and phenotypic data. Further details such as MRI modalities and scan parameters are available on the ABIDE website. This study was carried out in accordance with relevant guidelines and regulations.
Quality assessment and exclusion criteria before preprocessing. We used the following criteria to exclude subjects in ABIDE-I from this study. First, the subjects with missing anatomical or functional files were excluded. Second, all subjects from the imaging site Stanford were excluded as it is the only site with spiral image acquisition protocol. Third, all subjects from the imaging site Leuven-1 were excluded due to unknown repetition times for the functional scans. Fourth, to assess the quality of the raw images in ABIDE-I, we have used the information on raters' decisions available from the Preprocessed Connectome Project (PCP) 76 , and the subjects whose raw image quality was described as 'fail' by both the raters were excluded. Note that we did not exclude the subjects based on IQ or match the cohorts for IQ in order to ensure that the results of our analyses are generalizable to the typical ASD population 18,75 . After removing subjects based on the quality assessment (QA) checks Table 3. Agreement of results from node-level network analysis to TMS/tDCS literature. The list of the target brain regions that show improvement in clinical symptoms associated with ASD following TMS or tDCS procedure, the corresponding ROIs in the Schaefer atlas that show significant between-group differences in node-level network measures and the network measures that capture the differences (Forman-Ricci curvature (FRC), clustering coefficient (CC)). www.nature.com/scientificreports/ and exclusion criteria described above, we were left with 494 subjects in the ASD group and 520 subjects in the TD group (Supplementary Table S1).
Raw fMRI data preprocessing. We used the CONN functional connectivity toolbox 69 to process the rs-fMRI data from ABIDE-I. Figure 1 is a schematic summarizing the processing pipeline for rs-fMRI data used in this study. We have created a protocol video providing a visual guide to rs-fMRI preprocessing using CONN toolbox which is available at: https:// youtu. be/ ch7-dOA-Vlo. We performed motion correction, slice-timing correction, outlier detection, and structural and functional segmentation and normalization. First, the functional images were co-registered to the first scan of the first session. The SPM12 realign and unwarp procedure 77 was used to realign and motion correct the images using six rigid body transformation parameters: three translations in x, y and z directions, and three rotations namely pitch, yaw and roll. Second, the SPM12 slice-timing correction procedure 78 was used to temporally align the functional images. Third, Artifact Detection Tools (ART)-based outlier detection was performed where acquisitions with framewise displacement greater than 0.5 mm or global BOLD signal changes greater than 3 standard deviations were marked as outliers. Fourth, segmentation and normalization 79 was carried out to normalize the images into the standard Montreal Neurological Institute (MNI) space, and then, segment the brain into grey matter, white matter and cerebrospinal fluid (CSF) areas. Raw T1-weighted volume of the anatomical image and mean BOLD signal of the functional images were used as reference in this step. Subjects with bad image quality and signal dropouts in their scans or subjects with registration or normalization errors were excluded from further analysis.
After the spatial preprocessing of the raw rs-fMRI scans, the BOLD time-series associated with each voxel was extracted using the CONN toolbox. Next, we performed temporal preprocessing or denoising using the CONN toolbox to further reduce physiological or motion effects from the BOLD time-series. First, we implemented anatomical component-based noise correction procedure (aCompCor), to simultaneously remove 5 potential noise components 80 each from white matter and CSF areas, 12 potential noise components from estimated subject motion parameters and their associated first-order derivatives 81 , and 1 noise component from each of the identified outlier scans (scrubbing) 82 in a single linear regression step. Second, a high-pass filtering was performed to remove temporal frequencies below 0.008 Hz from the BOLD time-series.
Quality assessment and exclusion criteria after preprocessing. After preprocessing the raw fMRI data, we applied the following criteria to exclude participants from the analysis. Subjects were excluded if the FC distribution deviated significantly from normal distribution, or if the FC distribution showed noticeable distance dependence 83 . We additionally excluded subjects that showed a noticeable correlation between quality control (QC) variables and FC values, or if the QC-FC correlations showed a noticeable distance dependence 83 . After removing subjects based on these exclusion criteria, we were left with 395 subjects in the ASD group and 425 subjects in the TD group (Supplementary Table S1). The FC matrices of these remaining 820 subjects were used for network analysis. The demographic and clinical information for these subjects from ABIDE-I included in our study is summarized in Table 1.

Atlas-based definition of nodes and functional connectivity.
A widely-used approach for defining nodes in functional connectivity networks (FCNs) is to group closely related neighboring voxels into cortical parcels, in order to obtain nodes with interpretable neurobiological meaning 84 . Furthermore, the use of brain parcellations also reduces the computational load of further analyses. In this study, we used a predefined cortical parcellation atlas by Schaefer et al. 56 , which is based on a gradient-weighted Markov random field approach. While the Schaefer atlas is available at multiple resolutions, we considered the resolution that parcellates the brain into 200 distinct regions of interests (ROIs) wherein each hemisphere comprises 100 ROIs. In this parcellation, each ROI belongs to one of seven resting state networks (RSNs), namely, 'visual' , 'somatomotor' , 'dorsal attention' , 'salient ventral attention' , 'limbic' , 'control' , and 'default' . Using the CONN toolbox, the time series of each ROI was computed as the average of the time series of all the voxels that it contains. Subsequently, Pearson correlation coefficient between the time series of every pair of ROIs was calculated in the CONN toolbox, which resulted in a 200 × 200 FC matrix for each subject.

Construction of Sparsity-based functional connectivity networks.
In the preceding subsection, we described the FC matrix which is a correlation matrix that can be represented as a complete, weighted and undirected graph wherein the ROIs correspond to the nodes and the weights of edges are given by the correlation values between ROIs. The construction of the FCN from the FC matrix of a subject includes two steps, namely maximum spanning tree (MST) construction and sparsity-based thresholding. First, to extract the most important edges from the FC matrix, we constructed its MST using Kruskal's algorithm 85 . The MST is a spanning tree of the weighted graph with maximum edge weight. Note that the MST for a weighted graph with n nodes is an acyclic graph (more precisely, a tree) with ( n − 1 ) edges which is always connected. Second, we used sparsitybased thresholding, wherein edges are iteratively added to the MST in decreasing order of their correlation values, until a resulting network with the desired sparsity was obtained. Further, the resulting network with desired sparsity was binarized by ignoring the edge weights before proceeding to compute the network properties 26,86 .
Evidently, this choice of MST construction followed by sparsity-based thresholding to generate the FCNs ensures that the constructed networks for different subjects are connected and have the same number of edges. Such networks enable direct mathematical comparison of global and local network properties across subjects 26,87,88 . We remark that this choice of MST followed by sparsity-based thresholding to construct FCNs from rs-fMRI images has been used earlier by Achard et al. 70  www.nature.com/scientificreports/ As there is no rationale for using a specific graph density, previous studies 26,28,31 have explored network properties across a range of graph densities. In this work, we have studied the network properties over a wide range of graph densities between 0.02 or 2% edges and 0.5 or 50% edges, with an increment of 0.01 or 1% edges. Thus, for each of the 820 subjects from ABIDE-I considered in this study, we have constructed 49 unweighted and undirected networks. In other words, we have generated 820 × 49 FCNs for 820 subjects across 49 graph densities or thresholds for this study, and the constructed networks are made publically available in our GitHub repository.

Network based analysis and post-hoc analyses.
In this section we describe the methodology used to analyze and compare the resting-state FCNs of individuals with ASD and TD individuals constructed as mentioned in the preceding section. First, we performed global and node-level network analysis to compare the FCNs in the ASD group and the TD group. Second, we used Neurosynth-based meta-analysis decoding 89,90 to determine the behavioral relevance of the results of our node-level network analysis, and studied the relationship between node-level network measures and relevant scores of symptom severity in ASD. Third, we assessed the agreement of the results of our node-level network analysis against results reported in non-invasive brain stimulation (NIBS) studies with Transcranial Magnetic Stimulation (TMS) and transcranial Direct Current Stimulation (tDCS).
Global and node-level network analysis. As mentioned in preceding subsection, we constructed 49 unweighted and undirected networks with varying sparsity from the FC matrix corresponding to each subject, and thereafter, each of the 49 networks for a subject was characterized by computing discrete Ricci curvatures and other network properties. Specifically, we have focused here on two discrete Ricci curvatures, namely Forman-Ricci curvature (FRC) 39,41,42 and Ollivier-Ricci curvature (ORC) 40 . Notably, the two discrete Ricci curvatures are naturally defined for edges in a network and capture different aspects of the classical Ricci curvature 42 . Moreover, we have also explored here several standard global network measures including average clustering coefficient, modularity 91 , average shortest path length, average node betweenness centrality, global efficiency 92 and average local efficiency 92 . In Supplementary Information, we describe the different global and local network measures employed here to characterize the FCNs.
To compare the global properties of the FCNs across the two groups (ASD versus TD), we first computed the average FRC of edges, average ORC of edges and seven other global network measures (including average clustering coefficient, modularity, average shortest path length, average node betweenness centrality, global efficiency and average local efficiency), for each of the 820 × 49 networks corresponding to the FC matrices of 820 subjects across 49 graph densities. To compare the node-level properties of the FCNs across the two groups (ASD versus TD), we computed the node FRC and node ORC for each of the 200 nodes in each of the 820 × 49 networks corresponding to the FC matrices of 820 subjects across 49 graph densities. Note that the node Ricci curvature is defined as the sum of edge Ricci curvatures for the edges incident on that node 42 (see Supplementary  Information). Additionally, we computed two standard network measures, namely node clustering coefficient and node betweenness centrality.
The computer codes for FRC and ORC are made publically accessible via a GitHub repository. The other global network measures mentioned above for FCNs were computed using the Python package NetworkX 93 . Furthermore, the statistical tests were performed in Python packages SciPy 94 and statsmodels 95 .
Neurosynth meta-analysis decoding. We used Neurosynth meta-analysis decoding 89,90 to interpret the results of the node-based network comparisons, in terms of their behavioural relevance. Corresponding to each nodelevel network measure studied here, we identified a set of nodes (ROIs) that showed significant differences between the ASD and TD groups. For a set of nodes with significant between-group differences for a network measure, we used Neurosynth meta-analysis tool to find terms related to cognition, perception and behavior corresponding to the centroid coordinates of each ROI in the set. Further, we partitioned the set of identified ROIs which show significant between-group differences, by the 7 RSNs in the Schaefer atlas, and thereafter, the frequency counts of the terms associated with the subset of identified ROIs in a particular RSN were calculated and the statistical significance of these frequency counts was determined. This was done for each RSN separately, to identify those terms selectively associated with each of the 7 RSNs.
After determining the behavioral relevance of the brain regions with significant between-group differences in node-based network measures, we performed a post-hoc correlation analysis to measure the strength of the linear relationship between the values of the node-based network measure for each of the brain regions, and clinical scores related to symptom severity of the identified cognitive domains. We performed this analysis just for the ASD group. Specifically, we chose two clinical scores based on the Autism Diagnostic Interview-Revised (ADI-R) scoring 96 , namely ADI-R verbal and ADI-R social. We chose the ADI-R scores among all the possible clinical scores because they are available for the most number of participants ( n = 275 ) in the ASD group, and the ADI-R social and ADI-R verbal scores are appropriate means to capture symptom severity in autism compared to other clinical scores 97 .
Literature search for non-invasive brain stimulation studies in ASD. We performed a literature search to identify brain regions whose non-invasive stimulation were reported to result in improvement of ASDrelated symptoms. We first performed the literature search, to identify scientific papers reporting the effect of non-invasive brain stimulation (NIBS) on core symptoms of ASD, and then used results reported in these papers to identify those brain regions whose stimulation resulted in positive behavioral and cognitive outcomes. www.nature.com/scientificreports/ spectrum disorder' , ' Asperger's syndrome' , 'autism' and three major brain stimulation methodologies including 'Transcranial magnetic stimulation' , 'TMS' , 'transcranial direct current stimulation' , 'tDCS' , 'transcranial alternating current stimulation' , 'tACS' . The search was performed in October 2021 and the exact details regarding the search query are provided in Table 2. The PubMed search returned 235 articles. We followed a three-stage procedure to further refine the list of 235 articles returned by the PubMed search. First, we checked for papers missing from the corpus generated by the PubMed search, by scanning review articles on the use of NIBS methods to study ASD. We also searched these review articles for potential databases of NIBS experiments on ASD. Second, we filtered the articles based on title and abstract, based on relevance. We defined relevance according to the following criteria. The inclusion criteria were: (1) studies on ASD populations, (2) studies that have used NIBS techniques namely, TMS (and its variants such as rTMS), tDCS and tACS, (3) studies that have investigated the effect of NIBS on the core behavioral and cognitive symptoms of ASD, and (4) studies that are peer-reviewed. The exclusion criteria were: (1) review articles, (2) articles presented in languages other than English, (3) studies that did not perform NIBS, (4) studies that investigate new protocols for NIBS, (5) studies that report no positive effects in ASD symptoms post NIBS, (6) studies whose target areas for NIBS were not clearly reported, and (7) articles without access to full-text. Third, we classified the articles based on their stimulation technique (TMS/ tDCS/ tACS) and checked the full text of the articles for relevance, according to the same criteria as above. This process yielded 19 eligible articles for TMS, 12 eligible articles for tDCS and zero articles for tACS.
We identified Barahona-Corrêa et al. 98 as a database of TMS studies in ASD published before 2018, with data collection guided by preferred reporting items for systematic reviews and meta-analysis (PRISMA) 99 . Similarly, we identified García-González et al. 100 as a database of tDCS studies in ASD published before August 2019, also guided by PRISMA data collection. We utilized the data presented in these two databases along with the data that we extracted from the eligible articles in our corpus, such as author, publication year, DOI, number of participants, gender distribution, mean age, intellectual abilities, stimulation methodology and parameters, target areas, stimulation schedule, behavioral and cognitive outcome measures, behavioral and cognitive results, and any adverse reactions for the experiment group and the control group (if applicable). All the data collected are provided as Supplementary Tables S5 and S6. From these data, we identified the set of brain regions whose stimulation using these NIBS methods on individuals with ASD resulted in positive cognitive and behavioral outcomes.
Estimating overlap between regions identified in NIBS studies and node-level network analysis. We estimated the overlap between the sets of regions identified from literature search of NIBS studies and the sets of regions revealing ASD-related differences in node-level network measures. The target areas described in the NIBS studies were cortical regions in the brain that are specified by their respective Brodmann areas 68 while we identified node-level differences in areas of the Schaefer 200 atlas. We used the MRIcron tool 101 to map each of the Brodmann areas to Schaefer ROIs, by identifying the Brodmann area encompassing the MNI centroid coordinates of each Schaefer ROI 73 . The mapping from Schaefer ROIs to the Brodmann areas is presented in Supplementary Table S7. Next, we compiled the set of Brodmann areas that serve as target areas from the eligible NIBS experiments and have shown a positive outcome, either behavioral or cognitive, as a result of stimulating that region. We then identified the set of Schaefer ROIs that were mapped to these Brodmann areas. From this set of Schaefer ROIs, we found the subset that yielded significant ASD-related differences according to the graph Ricci curvatures namely, FRC and ORC, as well as for clustering coefficient and node betweenness centrality.
Quantification and statistical analysis. For the global measures, we evaluated the differences between the two groups across the 49 graph densities in the range 2-50% considered in this study by using a two-tailed two-sample t test. For the node-level measures, we first computed the area under the curve (AUC) for a given node measure across the 49 graph densities considered in this study 28,86 . Thereafter, we used a two-tailed twosample t test to evaluate the differences between the two groups via AUCs of the node measures for each of the 200 nodes in the network. Further, we measured the relationship between the values of the node-based network measure and the ADI-R scores by computing the partial correlations, with age and gender as covariates. For the Neurosynth meta-analysis decoding, to determine statistical significance of these frequency counts, we calculated the frequency counts of the same terms associated with an equal size set of randomly selected surrogate ROIs, and thereafter, the z-score for the frequency counts of each term associated with the subset of original ROIs was calculated. Subsequently, the z-scores were converted into p-values assuming a normal distribution.
After each of the above-mentioned tests or computations, we used a false discovery rate (FDR) correction 102 to correct for multiple comparisons and control the occurrence of false positives. Note that the alpha for these FDR corrections was set to 0.05.

Data and code availability
Functional connectivity matrices and networks generated in our study, and all original codes are deposited on GitHub and are publically available at https:// github. com/ asama llab/ Curva ture-FCN-ASD as of the date of publication.