The Resting-State Causal Human Connectome is Characterized by Hub Connectivity of Executive and Attentional Networks

We demonstrate a data-driven approach for calculating a “causal connectome” of directed connectivity from resting-state fMRI data using a greedy adjacency search and pairwise non-Gaussian edge orientations. We used this approach to construct n=442 causal connectomes. These connectomes were very sparse in comparison to typical Pearson correlation-based graphs (roughly 2.25% edge density) yet were fully connected in nearly all cases. Prominent highly connected hubs of the causal connectome were situated in attentional (dorsal attention) and executive (frontoparietal and cingulo-opercular) networks. These hub networks had distinctly different connectivity profiles: attentional networks shared incoming connections with sensory regions and outgoing connections with higher cognitive networks, while executive networks primarily connected to other higher cognitive networks and had a high degree of bidirected connectivity. Virtual lesion analyses accentuated these findings, demonstrating that attentional and executive hub networks are points of critical vulnerability in the human causal connectome. These data highlight the central role of attention and executive control networks in the human cortical connectome and set the stage for future applications of data-driven causal connectivity analysis in psychiatry.


Introduction
Brain network interactions give rise to information processing and cognition (Bressler, 1995;Bullmore and Sporns, 2009;Friston, 2002;McIntosh, 2000). Brain networks have a non-random topological organization (Bullmore and Sporns, 2009), including both segregated modules and a small number of highly connected nodes (Achard et al., 2006;Eguíluz et al., 2005;Sporns et al., 2007;van den Heuvel and Sporns, 2013), the "hubs" of the connectome. These hubs coordinate the transfer of large amounts of information through brain circuits (Mišić et al., 2015(Mišić et al., , 2014van den Heuvel et al., 2012) and play a critical role in coordinating communication between disparate brain regions (Cole et al., 2013;Sporns, 2013Sporns, , 2012van den Heuvel and Sporns, 2011). Here we present data on a "causal connectome" derived from resting-state neuroimaging data using a data-driven causal discovery method that first estimates directly connected brain regions without the false positives produced by typical connectivity methods (Reid et al., 2019), then additionally estimates the direction of those connections. We describe the characteristics of this causal connectome, with a specific focus on the central hubs of the resting-state causal connectome. While networks with directional information have several designations in the literature (effective, causal, directed), for consistency we will refer to these networks as causal networks throughout this manuscript.
Hub-like connectivity can be characterized using measures of centrality, a set of metrics that quantify the capacity of a node in a graph to influence (or be influenced by) other nodes. Over 100 measures of centrality have been proposed (Jalili et al., 2015), only a subset of which are commonly applied to brain network analysis (Rubinov and Sporns, 2010;van den Heuvel and Sporns, 2013). The most common is degree centrality (the number of connections attached to each node), which has a simple and intuitive interpretation, but alone provides limited information. For example, a node might have relatively few connections, but still act as an important bottleneck for communication between many other nodes in the network. Such a node could be thought of as a city highway. While it might have relatively few direct connections (entrances and exits), the highway still serves as the quickest path between many locations in the city. To capture this type of connectivity, measures such as betweenness centrality (how often a node lies on the shortest path between two other nodes; Freeman, 1977) have been developed.
Different metrics for centrality are often highly correlated with each other , but some centrality metrics may be more appropriate for certain types of networks than for others. For example, while degree centrality is frequently used for characterizing structural brain hubs (Crossley et al., 2014;Rubinov et al., 2015), degree in correlation-based fMRI networks is biased to identify nodes that are part of the largest (number of nodes/regions) resting-state networks (RSNs) of the brain van den Heuvel and Sporns, 2013). Thus, early degreebased analyses often identified high-degree nodes in the default mode network (Buckner et al., 2009;Cole et al., 2010;Power et al., 2013;Tomasi and Volkow, 2011;van den Heuvel and Sporns, 2013), one of the brain's largest resting-state subnetworks. Because of this confound, groups using Pearson correlationbased network analyses of fMRI data have considered alternative methods for identifying functional hubs, including metrics that consider the diversity of between-RSN connections such as participation coefficient (Bertolero et al., 2018;Grayson et al., 2014;Power et al., 2013;Reber et al., 2021), and coactivation/connectivity over multiple cognitive tasks (Cocuzza et al., 2020;Cole et al., 2013;Crossley et al., 2014Crossley et al., , 2013Ray et al., 2020). These more recent analyses reveal a set of brain hubs distributed broadly through frontal, parietal, and temporal cortices. Notably, defining what constitutes "high" centrality is dependent on the distribution of centralities exhibited by the network, and there is no evidence that one-size-fits-all cutoffs for categorizing network nodes as hubs can be meaningfully applied to brain networks. For example, when a cutoff for defining hubs as proposed in (Guimerà and Nunes Amaral, 2005) was applied to the brain network examined in Power et al. (2013), only a single local hub was identified. As such, the distinction between what constitutes a hub as opposed to a non-hub is essentially arbitrary . Our analysis of hub connectivity will thus focus on continuous measures of centrality, providing a quantitative comparative measure of which brain networks or nodes are more hub-like than others in the investigated causal connectome.
Critical to our current investigation, most resting-state fMRI connectivity studies use Pearson correlations to estimate undirected connectomes, providing information about connected brain regions (hereafter adjacencies), but not the direction of these connections (hereafter orientations). There is a recognized need for network modeling methods that can accurately estimate adjacencies without the false positives inherent to correlation-based approaches (Reid et al., 2019;Smith et al., 2011), as well as estimating the direction or orientation of these edges (referred to as causal or effective connectivity (Ramsey et al., 2010;Reid et al., 2019;Smith, 2012). However, current methods for fMRI resting-state causal connectivity are limited (Ramsey et al., , 2010Sanchez-Romero et al., 2019;Smith et al., 2011). Granger causality (Granger, 1969) attempts to recover causal influences using time-lagged regressions. Smith et al. (2011) found that several variations of Granger causality have negligible accuracy in detecting adjacencies or orientations in simulated fMRI data, and Sanchez-Romero et al.
(2019) tested two additional variations of Granger causality (multivariate Granger causality; Barnett and Seth, 2014; autoregressive modeling with permutation testing; Gilson et al., 2017), finding that these more recent methods also had low precision for adjacencies and orientations in the presence of realistic noise. GIMME (Gates and Molenaar, 2012), a group-level algorithm that also uses time lags to infer causality, achieves better performance, but is computationally intensive and can only scale to a small number of brain regions . Dynamic causal modeling (DCM; (Friston et al., 2003) was originally designed for task-based fMRI data, but the stochastic DCM (Li et al., 2011) and spectral DCM (Friston et al., 2014) variants can be applied to resting-state fMRI data (albeit only for a small number of brain regions due to computational demands). Frässle et al. (2021) recently proposed a highly scalable variant of regression DCM (Frässle et al., 2017) for resting-state fMRI data that is scalable enough to support whole-cortex analyses. While the connectomes generated by this newly proposed method appear to have face validity (Frässle et al., 2021), a quantitative analysis of large-scale restingstate connectivity patterns using regression-DCM has not yet been completed. Instead, the current report focuses on a distinct set of methods grounded on Bayes networks, which have promise for uncovering causal connectivity from fMRI (Mumford and Ramsey, 2014;Ramsey et al., 2014;, and for dealing with the high dimensionality of whole-cortex data (Ramsey et al., 2017).
Unlike GIMME, spectral DCM, and stochastic DCM, Bayes net methods are highly scalable, and unlike regression DCM, these methods do not require specification of priors or hemodynamic response function and make no assumptions about the physiology giving rise to the observed hemodynamic signal. Suitable combinations of causal discovery-based methods can achieve near-perfect precision and recall in simulated fMRI data (Hyvärinen and Smith, 2013;Ramsey et al., 2014), even for networks with feedback cycles ).
In the current study, we capitalize on these recent advances in causal discovery machine learning to build whole-cortex causal connectomes from single-subject resting-state fMRI data. We apply a variation of a previously proposed Sanchez-Romero et al., 2019) two-step causal discovery framework that breaks the connectome computation into separate adjacency and orientation steps. For convenience we refer to our approach as GANGO (Greedy Adjacencies and Non-Gaussian Orientations). GANGO first estimates whole-cortex adjacencies using Fast Greedy Equivalence Search (FGES; Ramsey et al., 2017). FGES is a parallelized version of Greedy Equivalence Search (Chickering, 2002), a Bayes Network method with high sensitivity for detecting adjacencies, but poor accuracy for orientations in simulated fMRI data (Smith et al., 2011). Thus, we follow this initial adjacency search with a pairwise edge orientation algorithm that exploits non-Gaussian information in the hemodynamic signal (Hyvärinen and Smith, 2013), shown in Sanchez-Romero et al., 2019) to have high precision and recall for determining edge orientation of the Smith and colleagues (2011) simulations. We thus obtain, on a single-subject basis, a whole-cortex graph summarizing dominant causal connectivity between brain regions. We focus our investigation on the hub structure of this novel causal connectome. Overall, we demonstrate hub-like causal connectivity profiles of the dorsal attention network, frontoparietal network, and cingulo-opercular network.

Resting-State fMRI Acquisition and Preprocessing
Structural (T1 and T2 images, required for preprocessing functional neuroimaging data) and functional MRI data were collected at Washington University on the Siemens 3T Connectome Skyra scanner. Full details of the acquisition parameters for the HCP data are described in .
Each subject's resting-state data was collected over two days in four sessions (14:33/session; 1200 samples/session). In this study we analyzed only one day of data (two runs, individually z-scored and concatenated) to limit potential state influences on fMRI measures. Structural and functional data preprocessing is described in , and used version 3.21 of the HCP preprocessing pipeline. Structural data preprocessing consisted of bias field and gradient distortion correction, coregistration of T1/T2 images, and registration to MNI space. Cortical surface meshes were constructed using FreeSurfer, transformed to MNI space, registered to individual surfaces, and downsampled.
Functional MRI preprocessing consisted of gradient distortion correction, motion correction, EPI distortion correction, followed by T1 registration. Transforms were concatenated and run in a single nonlinear resampling to MNI space followed by intensity normalization. Data were masked by the FreeSurfer brain mask, and volumetric data were mapped to a combined cortical surface vertex and subcortical voxel space ("grayordinates") using a multimodal surface registration algorithm (Robinson et al., 2014) and smoothed with a 2mm FWHM Gaussian in surface space (thus avoiding smoothing over gyral banks). fMRI data were conservatively high-pass filtered with FWHM = 2000 s and cleaned of artifacts using ICA-FIX Salimi-Khorshidi et al., 2014). This filter was implemented as a weighted linear function in FSL v6.0.2, which was shown in  to not introduce Gaussian trends into the data (unlike e.g., Butterworth filters or the built-in SPM filter).
Artifact components and 24 motion parameters were regressed out of the functional data in a single step, producing the final ICA-FIX denoised version of the data in CIFTI ("grayordinates") space (Glasser et al., 2016b) that was used in subsequent analyses.

Construction of Whole-Brain Causal Connectomes
Our analysis pipeline began with n = 442 sets of fully preprocessed, multi-modally surface registered, ICA-FIX denoised fMRI data provided by the HCP consortium. We parcellated cortex surface vertices into 360 regions using a recently developed multimodal parcellation (Glasser et al., 2016a).
We implement a computational strategy to define causal connectomes on a per-subject basis using a two-step process we refer to as GANGO (Greedy Adjacencies and Non-Gaussian Orientations).
This approach is motivated by previous work Sanchez-Romero et al., 2019;Smith et al., 2011) indicating that 1) Bayes net algorithms such as PC (Spirtes et al., 2001) and Greedy Equivalence Search (Chickering, 2002) provide a highly precise solution to identify nodal adjacencies (but not orientations) in simulated fMRI data, and 2) pairwise orientation algorithms based on data skewness can accurately identify edge orientations in simulated fMRI data. In the first step, the GANGO approach defines nodal adjacencies (connected regions) using Fast Greedy Equivalence Search (FGES; Ramsey et al., 2017), a parallelized and highly scalable version of GES. This algorithm finds a sparse set of directed and undirected connections between continuous variables by minimizing a penalized likelihood score over the entire graph, typically scored using the Bayesian Information Criterion (BIC; Schwarz, 1978). FGES proceeds in two stages, first greedily adding edges until the BIC stops improving, then greedily removing edges until the BIC stops improving (Ramsey et al., 2017). While FGES has not commonly been used for analysis of empirical neuroimaging data, this method was applied (in combination with direct stimulation) to test causal connectivity patterns of the amygdala (Dubois et al., 2020) with promising initial findings for mapping human emotion networks. We computed FGES with causal-cmd v1.2.0 (https://bd2kccd.github.io/docs/causal-cmd/) using default parameters (BIC penalized likelihood score, penalty discount = 1 corresponding to the classic BIC score). GES has been shown in simulations to obtain highly accurate estimates of nodal adjacencies, but relatively inaccurate orientations (Smith et al., 2011). Therefore, we made the FGES-derived graph undirected by symmetrizing it across the diagonal.
The GANGO approach then orients these undirected edges using an estimate of the direction of causal effect based on pairwise likelihood ratios under the linear non-Gaussian acyclic model (Hyvärinen and Smith, 2013). Hyvärinen and Smith (2013) provide an explanation of how non-Gaussian information can be used to orient edges between pairs of variables. Assuming x -> y, both variables will have large values in cases where x is large (but x will not necessarily take on large values when y is large). Due to regression towards the mean, the value of x must typically be larger than that of y. Cumulant-based approaches such as RSkew, the method used here, calculate pairwise contrasts that magnify extreme values of either x or y, allowing determination of the most likely causal direction. For an in-depth explanation, we refer the reader to (Hyvärinen and Smith, 2013). We used the authors' MATLAB implementation of RSkew (https://www.cs.helsinki.fi/u/ahyvarin/code/pwcausal/; RSkew is method 4), an outlier-robust skew-based measure of the pairwise likelihood ratio (Hyvärinen and Smith, 2013), which has shown to generate optimal estimates of causal direction in simulated fMRI data Sanchez-Romero et al., 2019).
Since the RSkew orientation method requires that the data be non-Gaussian, and specifically skewed, we first tested whether the resting-state data met these assumptions. In the first analysis, we tested the data for non-Gaussianity using an approach from Ramsey et al. (2014). We calculated Anderson-Darling tests (MATLAB adtest function) over each participant's distribution of fMRI data and returned a test statistic and corresponding p-value for the alternate hypothesis that the data are non-Gaussian. In the second analysis, we tested whether the fMRI data were significantly skewed with reference to Gaussian data using an approach from Sanchez-Romero et al. (2019). We calculated the mean and standard deviation of every participant's distribution of fMRI data, and simulated Gaussian data with the same mean and SD for use as surrogate data (100 permutations per subject). For every run (out of 100), we calculated the skewness of the surrogate Gaussian data. For each subject, we then used the 95th percentile of this surrogate distribution as a critical value. We then compared the skewness of the participant's distribution of fMRI data to the critical skewness derived from Gaussian surrogate data.
Since both positive and negative skewness values indicate skewed data, for all analyses we used the absolute value of skewness. Figure 1. Summary of the strategy we employed to build single-subject (n = 442) causal connectivity graphs for further analysis, and the analyses we ran to characterize the hub structure of these RSNs.

Resting-State Network Connectivity Statistics
For each subject, we categorized brain regions into 12 resting-state networks ( RSNs (the total number of out-of-RSN connections). Then, by dividing that total by 11, we arrive at a suitable null hypothesis of equal inter-RSN connectivity (i.e., that the RSN shared out-of-RSN connections equally between the other 11 RSNs). For each of the other 11 RSNs, we tested the actual number of shared between-RSN connections against the null hypothesis of equal connectivity to all 11 RSNs using Wilcoxon signed-rank tests, which we corrected for multiple corrections using the false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995), thus establishing whether pairs of RSNs were significantly connected. To clarify the direction of causal connectivity between significantly connected pairs of RSNs, we calculated the proportion of causal connections from each RSN to each other RSN, then compared that proportion against a null hypothesis of 50% (i.e., that the connections between the pair of RSNs A & B are equally from A to B and from B to A) using Wilcoxon signed-rank tests (FDR-corrected for multiple comparisons). For each pair of RSNs found to be significant with reference to a null hypothesis of equal (i.e., random) inter-RSN connectivity, we calculated the effect size of the test for significant shared connections using Cohen's d, and masked inter-RSN connections that did not meet at least the threshold for a small effect size (Cohen's d >= 0.2). We plot this result in the form of a force-directed graph, by first drawing connections between pairs of significantly connected networks. If we were able to determine a significant direction of connectivity, this connection was unidirectional; otherwise, the connection was drawn bidirectional. For significantly one-directional connections, we also recorded the proportion of connections going in the significant direction. This analysis clarified a) which RSNs are significantly connected with at least a small effect size, and b) which RSNs significantly send (vs. receive) information to (vs. from) other RSNs. We supplemented this with an analysis of participation coefficient (a measure of the diversity of RSNs a node connects to), calculated separately for incoming and outgoing connections. Participation coefficient Pi of node i is calculated according to the equation , where  is the number of connections between node i and module s and Ki is the total degree of node i (Guimerà and Nunes Amaral, 2005). For nodes that connect entirely within their own module, this equation results in a participation coefficient of zero, and for nodes that connect homogeneously over all modules the participation coefficient will approach one (Guimerà and Nunes Amaral, 2005;Power et al., 2013). We computed participation coefficients using the part_coeff function in the Brain Connectivity Toolbox (Rubinov and Sporns, 2010). We compared nodal participation coefficients in 12 established RSNs (Ji et al., 2019) that have been validated for our current parcellation using Friedman tests (since network is a within-subject measure) and conducted post-hoc pairwise comparisons with control for multiple comparisons using the Nemenyi test.

Centrality Distributions in Human Cortex
For each causal graph, we calculated nodal (n = 360) centrality statistics using indegree (number of incoming connections), outdegree (number of outgoing connections), and betweenness centrality (number of shortest paths the node participates in) (Figure 1c). To statistically quantify whether centrality-based cortical hubs existed based on significantly heavy-tailed centrality distributions, we generated a reference set of 1000 random directed graphs with the same number of nodes (360) and connections as the cortical causal graphs (Figure 1c). Specifically, for each run (of 1000), we chose the exact number of connections by drawing a random number from a normal distribution with the same mean and standard deviation as the number of connections across subjects (mean n connections = 1452, SD = 107); thus, the surrogate graphs approximated the distribution of connection counts present in the actual data. Random graphs were created using the makerandCIJ_dir function from the Brain Connectivity Toolbox (Rubinov and Sporns, 2010), which creates a random directed (causal) graph with a specified number of connections. We then used Wilcoxon rank-sum tests to compare the skewness of centrality distributions (indegree, outdegree, betweenness) of the causal connectivity graphs to the random graphs. We additionally applied this permutation analysis to categorize cortical nodes as hubs or non-hubs, based on whether a node's centrality exceeded the 95 th percentile of the surrogate distribution (i.e., whether a node was significantly in the highly central tail of the distribution). These binary masks of regions thus categorized as significant hubs are provided in the supplement.

Resting-State Network Differences in Centrality-Based Hubs
We compared the average nodal centralities (degree, betweenness) in 12 established RSNs (Ji et al., 2019) that have been validated for our current parcellation (Figure 1d), to obtain a continuous ranking of which networks are the most "hub-like" in the causal connectome. For each subject, we calculated the average nodal indegree, outdegree, total degree, and betweenness centrality for nodes within each of these 12 RSNs using the MATLAB centrality function. We compared centrality across the 12 RSNs using Friedman tests and conducted post-hoc pairwise comparisons with control for multiple comparisons using Nemenyi tests. Additionally, we compared indegree and outdegree within each of the 12 RSNs using Wilcoxon signed-rank tests, and FDR-corrected the resulting p-values for n = 12 multiple comparisons.

Network Vulnerability to Targeted and Random Attack
To clarify the functional role of hubs in the causal cortical network, we subjected each subject's causal connectome to a targeted attack analysis. Broadly, virtual attack analyses proceed by deleting nodes from the network and recording some measure of whole-network fitness as nodes are removed. The dependent value in this analysis is chosen to be a measure of network functional integration, to characterize nodes that are most critical for network communication (Rubinov and Sporns, 2010). As discussed in (Rubinov and Sporns, 2010), the most common measure of network functional integration is the characteristic path length (Watts and Strogatz, 1998). However, characteristic path length cannot be meaningfully computed on networks with disconnected nodes, as disconnected nodes are defined to have infinite path length. Thus, authors have argued that global efficiency (the inverse shortest path lengths in the network, defined to be zero for disconnected nodes (Latora and Marchiori, 2001)) is a superior measure of network functional integration (Achard and Bullmore, 2007;Bassett and Bullmore, 2006;Rubinov and Sporns, 2010). As such, network global efficiency is typically used as the dependent value in published virtual attack analyses (Crossley et al., 2014;Lin et al., 2018;Lo et al., 2015;van den Heuvel et al., 2018;van den Heuvel and Sporns, 2011).
To assess cortical vulnerability at the RSN level, we deleted the nodes of each RSN (one at a time) from the individual subject-level cortical causal connectivity graphs, and we recorded changes in connectome communication efficiency as percent change in global efficiency (Figure 1e), calculated using the efficiency_bin function in the Brain Connectivity Toolbox (Rubinov and Sporns, 2010). We supplemented this targeted attack analysis with a random attack analysis, by deleting randomly selected nodes rather than specifically targeting nodes from one RSN at a time. To minimize order effects in the virtual lesion analysis we ran the analysis 10 times for every subject, deleting the nodes within each RSN in a random order each time, then we took the mean of the efficiency loss over the 10 runs. This resulted in a set of 13 loss-of-efficiency curves per subject that quantified how strongly communication was impaired as successive nodes from each RSN were deleted. We took the pointwise derivative of each subject's loss-of-efficiency curves (for each of 13 deletion schedules) and compared the average pointwise slope between RSNs (plus random deletion) to quantify which RSN resulted in the most rapid loss-ofefficiency. RSN deletion slopes were statistically compared using a Friedman test with post-hoc significance testing using the Nemenyi test. Additionally, we examined which nodes had the greatest overall impact on connectome loss-of-efficiency by deleting (in single subjects) each node (of 360) one at a time and recording the change in overall global efficiency. Results

Individual Subjects Functional Neuroimaging Data are non-Gaussian and Skewed
Anderson-Darling tests applied to every subject's concatenated fMRI time series found significance in 438/442 subjects (99%), rejecting the null hypothesis that the data are Gaussian. In a more stringent permutation-based analysis, we compared the skewness of each participant's fMRI time series to the skewness of random Gaussian data. This analysis revealed that 421/442 (>95%) of participant's fMRI time series were significantly (p < .05) more skewed than surrogate Gaussian data. At the group level, data were significantly more skewed than random Gaussian data as well (p < .001). Thus, we conclude that the functional neuroimaging data analyzed in the current report are non-Gaussian, specifically skewed, and therefore it is appropriate to apply the skew-based orientation step to the analysis.

Whole-Cortex Causal Graphs are Sparse but Well-Connected
The resulting whole-cortex causal connectomes were sparse, containing only ~2.25% of all possible connections (360 parcels: mean n connections = 1452). Nevertheless, the graphs were wellconnectedin most graphs (93.7%), every node was connected to at least one other node, and of the 6.3% of graphs that contained disconnected nodes, each graph had very few disconnected nodes (median = 1, max = 2). Thus, despite the sparsity of the cortical causal connectivity graphs, these graphs appear to capture global causal patterns of connectivity.

Diversity of Inter-Network Connections Highlights Hub Roles of Multiple Brain Networks
To summarize the overall connectivity structure of our cortical causal graphs, we categorized each of 360 cortical nodes into 12 resting-state networks (RSNs) (Ji et al., 2019). For a visual of these 12 RSNs plotted on the cortical surface, see Figure 1A. We first examined patterns of connectivity between these large-scale brain RSNs. We established whether each RSN shared a statistically significant proportion of connections with each other RSN. We additionally tested for a preferred direction of connectivity between significantly connected pairs of RSNs.
The result of this analysis is plotted as a force-directed graph, with connections drawn between RSNs that were found to be significantly connected. Connections are shown as unidirectional if rank-sum testing indicated a preferred direction of connectivity, and connections are shown as bidirectional if ranksum testing was unable to determine a preferred direction of connections (Figure 2a). A clear hubperiphery structure was apparent. We found that visual RSNs 1 and 2 were bidirectionally interconnected, and that visual RSNs projected to the dorsal attention network. The dorsal attention network projected to multimodal association networks (posterior and ventral) and the frontoparietal network. The frontoparietal network was situated in the center of the graph, being the most highly connected RSN and sending information to the ventral attention/language, limbic/orbito-affective, and default mode RSNs, while bidirectionally sharing connections with the cingulo-opercular network. This network connectivity diagram also demonstrated significant bidirectional connectivity between auditory network and the ventral attention/language network, as well as showing strong connectivity between somatomotor and auditory network, which is to be expected since these networks are spatially adjacent and are strongly interconnected (so strongly that they are frequently merged into the same network in the literature; Ji et al., 2019).
We supplemented this analysis by an examination of the average participation coefficient within each RSN. Since we used causal connectivity graphs, we calculated participation coefficient separately using outgoing and incoming connections, resulting in two summary measures per RSN, per subject (hereafter out-part and in-part, respectively). We found that across subjects, the twelve RSNs differed

The Causal Connectome Has a Heavy-Tailed Centrality Structure
Thus far we have clarified patterns of inter-RSN connectivity in the causal connectome. From here, we examined the most important hubs of the cortical causal connectome using two common centrality metrics: degree centrality, which characterizes highly connected nodes, and betweenness centrality, which characterizes nodes that lie on many shortest paths between other nodes and thus facilitate efficient network communication (Rubinov and Sporns, 2010). If nodes are to be considered hubs based on centrality metrics, the distribution of those metrics should be heavy tailed, with a minority of highly central nodes. We found that the centrality distributions of the obtained causal graphs were indeed heavy-tailed, with most nodes having very low centrality and relatively few nodes having very high causal centrality (Figure 3a,b,c,d). Centrality values of the causal connectivity graphs were significantly more heavy-tailed than equally connected random comparison graphs, as confirmed by comparing the skewness of the indegree, outdegree, total degree, and betweenness distributions of these two sets of graphs (rank-sum tests, all p < .001). We additionally leveraged this analysis to produce categorical labels of individual brain regions as hubs if their centrality exceeded the 95 th percentile of the surrogate distribution. These binary labels almost exclusively designated regions of the lateral and superior parietal cortex as hubs, along with some frontal nodes (Supplement).

Networks
The 12 RSNs differed in average indegree, outdegree, and total degree (Friedman tests, all p < .001). Most relevant for the current work, post hoc testing for RSN degree differences (Nemenyi tests) demonstrated that the dorsal attention and frontoparietal networks had significantly higher indegrees, outdegrees, and total degrees than the other 10 RSNs (all p < .001). Similarly, comparison of average betweenness centrality across the 12 RSNs (Figure 3j) established that frontoparietal nodes participated in the greatest number of efficient paths, followed by dorsal attention nodes; these RSNs had higher betweenness centrality averages than the other 10 RSNs, and frontoparietal had significantly higher betweenness than dorsal attention (all p < .001). The cingulo-opercular network had the third highest average betweenness centrality scores, despite having only modest degree centrality, thus suggesting that while cingulo-opercular regions might not be the most highly connected regions in cortex, these regions are nevertheless particularly important for cortical communication. Note that while Power and colleagues (2013) showed that degree-based hubs in Pearson correlation networks are confounded by the size of the functional communities the nodes belong to (i.e, the number of nodes in each RSN), in a critical control analysis we did not find significant correlations between community size and degree or betweenness centrality (Figure 3k,l,m,n), suggesting that our measures of causal centrality cannot be ascribed to the size of the RSNs in our analysis.

Executive and Attentional Networks Equally Send and Receive Connections
Since causal graphs separate incoming and outgoing causal connections, we were additionally able to assess whether each RSN primarily sent or received information. Most RSNs could be characterized as either primarily "senders" (visual, somatomotor), or as primarily "receivers" (cinguloopercular, auditory, default mode, posterior/ventral multimodal, and limbic/orbito-affective). However, a small number of RSNs were found to send and receive equal numbers of connections (frontoparietal, dorsal attention, ventral attention/language; Figure 3i).

Figure 3. The Most Central Hubs of the Causal Connectome Cluster in Executive and Attentional
Networks.

A,B,C,D: Median nodal centrality across n = 442 subjects, for each of 360 cortical nodes. Blue
histogram indicates the distribution of median centralities for 1000 equally connected random graphs.

K,L,M,N: Control analyses ruled out the possibility that cortical hubs could be explained by the number of parcels in the RSN that each node belongs to. Blue line indicates a least-squares regression fit, and blue
shading indicates the 95% confidence interval of the regression.

Executive and Attentional Hubs are Points of Causal Connectome Vulnerability
Our analyses up to this point have demonstrated that the frontoparietal, dorsal attention, and cingulo-opercular RSNs are highly central hubs in the causal connectome. Based on previous reports, it is likely that the identified central hubs of the causal human connectome are also points of system-level vulnerability to insult. To test this hypothesis, we conducted a series of simulated attacks on the causal connectomes presented in this study. For each RSN, we sequentially deleted nodes in that RSN from each subject's cortical graph, and measured loss-of-function via percent change in global efficiency (Latora and Marchiori, 2001;Rubinov and Sporns, 2010). Figure 4a shows the average of these network-level loss-of-function curves for all 442 subjects. As a summary measure of the impact of nodal targeted attacks on each RSN, we took the average pointwise derivative of the global efficiency loss curve for each subject and RSN (plus random deletion, as a control analysis; Figure 4b). Results indicated that the RSN loss functions differed significantly in average pointwise slope (Friedman test, p < .001). Post-hoc multiple testing (Nemenyi test) indicated that the frontoparietal network had the steepest loss-of-

C: Change in network global efficiency following deletion of individual cortical nodes. Color indicates
the loss-of-efficiency when each single node is deleted.

Comparison of GANGO Causal Connectivity Graphs with Pearson Correlation Graphs
To examine how the cortical causal human connectome compares to more typical connectivity analyses (Pearson correlation graphs), we ran the presented analyses using two sets of binarized correlation graphs. The first set was proportionally thresholded at a 15% cost (that is, each graph retained the 15% largest positive values). Proportional thresholding was chosen to improve stability of measures over absolute thresholds (Garrison et al., 2015), and the chosen 15% cost is in the middle of an ideal cost range for producing small-world graphs in Pearson correlation networks (Achard and Bullmore, 2007;Bullmore and Bassett, 2011) and as such is among the most typical thresholding procedures in the literature. The second set was thresholded to retain the same number of connections as the subjectspecific causal graph, thus matching the density of the causal connectomes exactly. Detailed results of these comparison analyses are presented in the Supplement.
Overall, this comparison suggests that Pearson correlation graphs emphasize the importance of sensory regions and motor cortex as cortical hubs, while causal connectivity graphs instead emphasize higher cognitive regions, in particular frontoparietal and cingulo-opercular networks (which did not exhibit hub-like connectivity in any analysis for Pearson correlation graphs). We also found that the sparser set of correlation graphs were very poorly connectedon average, over 40% of nodes were completely disconnected at this threshold. At this 2.25% (average) density, we also found that graphs were significantly less hub-like than surrogate random graphs, suggesting that centrality-based hubs break down at this level of sparsity in Pearson correlation graphs.
Furthermore, in 15% density Pearson correlation graphs we found that both degree and betweenness centrality were highly confounded by the size of the RSN nodes belonged to, unlike in the causal graphs. This was not the case for the sparser correlation graphs. These differences between Pearson correlation and causal connectomes were somewhat attenuated when using participation coefficient as a measure. However, in the sparse 2.25% density correlation graphs, we were unable to estimate participation coefficient for the two least-connected networks (ventral multimodal, limbic/orbitoaffective) due to extremely low levels of connectivity (zero in nearly all subjects). We also observed that applying the virtual lesion analysis to the 15% density Pearson correlation graph resulted in many RSNs increasing global efficiency when deleted. This effect was magnified for the sparser 2.25% density correlation graphs. This is likely due to the much greater incidence of completely unconnected nodes in thresholded Pearson correlation graphs, as opposed to the well-connected causal graphs we generated using the GANGO method. Note that unconnected edges have this effect because the global efficiency of a node is calculated as the inverse shortest path (number of edges) from that node to each other node in the network (Latora and Marchiori, 2001). The global efficiency of the network is then the average of all nodal global efficiencies. An unconnected node is defined to have an infinite path length and a global efficiency of zero (Rubinov and Sporns, 2010). Removing this zero from being averaged into the global efficiency score thus improves the network efficiency.

Discussion
Functional connectivity analyses of the human cortex typically use undirected connectivity estimates, derived by computing Pearson correlation coefficients between the time series of the hemodynamic signals of individual brain regions (nodes). The need to extend functional connectivity analyses to causal connectivity is recognized (Reid et al., 2019;Smith, 2012), but data-driven methods for calculating high-dimensional causal graphs within single subjects have not been properly tested yet. Here, we present an examination of the causal connectivity patterns of the human cortex using a two-stage causal discovery machine learning approach. Two-stage causal discovery methods break the graph creation process into separate adjacency search and orientation phases. Ramsey et al. (2014) (Chickering, 2002), which was shown to produce high precision for adjacency search in Smith et al. (2011).
The causal connectomes produced by the GANGO approach were quite different from those produced by the more typical Pearson correlations -despite very low density, these graphs were fully connected in nearly all cases. Furthermore, GANGO networks did not exhibit any relationship between RSN size (number of nodes) and degree centrality, unlike standard Pearson correlation-based graphs. This dependency arises in correlation-based graphs whenever the graph exhibits a modular community structure, as explained in Power et al. (2013). Specifically, in cases where graphs exhibit community structure, a node in a larger community would by definition have a higher chance to form more connections simply because connections are more common within communities than between them.
However, many of these connections will actually be indirect, and nodal correlations within a community will be high due to indirect connections. Bayes net methods, on the other hand, enforce sparsity wherever possible via a Markovian screening-off property and retain only direct connections while eliminating indirect connections (Spirtes et al., 2001). Thus, our results provide support for the viability of the GANGO approach for providing unbiased centrality measures from resting-state fMRI data.
Prominent hubs of the causal connectome overlap many regions previously identified by restingstate fMRI (Achard et al., 2006;Buckner et al., 2009;Tomasi and Volkow, 2011;van den Heuvel and Sporns, 2013;Zuo et al., 2012), with the GANGO method reliably recovering these network properties when applied on a single subject level. Importantly, control analyses indicated that nodal hub metrics (degree, betweenness centrality) were unconfounded by the size of the RSN that nodes belonged to. In contrast, degree and betweenness centrality were strongly confounded by RSN size for thresholded correlation graphs (Supplement). Overall, we found prominent causal connectivity hubs and points of vulnerability of the causal connectome in dorsal attention network (DAN), frontoparietal network (FPN) and cingulo-opercular network (COP), with each of these hub networks showing distinctly different connectivity profiles.
The dorsal attention network (DAN) exhibited theoretically interesting properties that contributed to its high level of connectivity. In our analysis, DAN had among the greatest diversity of connections with other RSNs (measured using participation coefficient), as well as having overall high connectivity (measured using degree centrality) and participating in many efficient paths (measured using betweenness centrality). Our analysis of the inter-RSN connectivity structure of the cortical causal connectome revealed that DAN owed its high centrality to its role in receiving information from visual networks, processing that information, and then transmitting information to multimodal association networks (posterior/ventral multimodal association networks, FPN). This is in line with long-standing evidence that DAN plays an important role in top-down visual selective attention (Corbetta and Shulman, 2002;Vossel et al., 2014). While we found that causal connections usually progressed from visual to dorsal attention networks, a large proportion of these connections still progressed in the opposite direction as well. Thus, our results support a role of DAN in top-down control over visual systems as well, providing further evidence that the dorsal attention network supports both bottom-up sensory integration and top-down attentional control (Long and Kuhl, 2018).
The cingulo-opercular network (COP) was found to mediate many efficient paths in the cortex (betweenness) and shared a large diversity of inter-RSN connections (participation coefficient) but did not have particularly high connectivity (degree). COP has a role in maintaining task sets, initiating goaldirected behaviors, and consolidating motor programs (Dosenbach et al., 2008(Dosenbach et al., , 2006Fair et al., 2007;Newbold et al., 2021). In our consensus network structure, we found that COP received significant connections from the somatomotor network, sent significant connections to the orbito-affective (reward) network, and was bidirectionally connected to FPN. The uncovered functional connectivity between the reward networks and COP is in line with evidence that COP has a role in coordinating the response of brain reward-related regions (Huckins et al., 2019), and the connectivity between COP and somatomotor networks corroborates evidence that COP plays a role in consolidation, planning, and plasticity of motor regions (Newbold et al., 2021). Connectivity between COP and somatomotor networks also increases through development and is linked to the development of improved cognitive control (Marek et al., 2015).
Finally, COP was found to be tightly bidirectionally connected with frontoparietal network, echoing evidence that these RSNs work together as dual cognitive control networks (Dosenbach et al., 2008;Fair et al., 2007;Gratton et al., 2018).
Across all analyses, we found a critical role of frontoparietal executive network (FPN) connectivity. This represents an important point of agreement between our causal network results and recent advances in understanding the role of the FPN in overall function, as well as an important result that we did not find with traditional Pearson correlation graphs (Supplement). Due to its central position, FPN shared the greatest diversity of inter-RSN connectivity in the consensus graph, including significant received connections from DAN, significant sent connections to the ventral attention, orbito-affective, and default mode networks, and bidirectional connections with COP. Previous studies demonstrate that FPN flexibly shifts its connectivity patterns to fulfill task demands, while still retaining high correlations with its resting-state connectivity (Cocuzza et al., 2020;Cole et al., 2013;Crittenden et al., 2016), and another recent study demonstrated that resting-state network connectivity in FPN predicts transfer of taskrelevant information through distributed brain circuitry (Ito et al., 2017). Overall, our finding that FPN nodes are, on average, the most highly central nodes in the causal human connectome is consistent with a theoretical role of FPN as a flexible executive coordinator of overall brain function (Assem et al., 2020;Dosenbach et al., 2006;Duncan, 2010;Fair et al., 2007;Marek and Dosenbach, 2018).
This is also consistent with recent control systems perspectives on brain connectivity, which have suggested that FPN has a particular role in shifting brain network configuration into difficult-to-reach cognitive states (Gu et al., 2015). The causal connectivity from FPN to default mode network (DMN) deserves particular attention. These RSNs are typically anticorrelated during cognition (FPN activity increases and DMN activity decreases), a finding that replicates in experimental paradigms including attention (Hellyer et al., 2014), working memory (Kelly et al., 2008;Murphy et al., 2020), and cognitive reasoning (Hearne et al., 2015). We found that the connectivity between FPN and DMN preferentially flows from FPN to default mode network, suggesting that FPN might have an executive role in "turning off" or inhibiting the DMN in response to the need for control. The revealed causal hub role of the FPN is among the most important contributions of this study, as typically thresholded Pearson correlation graphs do not show hub-like connectivity in FPN (Supplement), despite ample theoretical and experimental evidence that the frontoparietal network is critical for overall organization and control of the connectome.
The current investigation delivers a powerful new framework for quickly computing (~30 s per connectome on a personal laptop with six cores and 32 Gb Ram) high-dimensional causal connectivity graphs from observed brain data as well as providing important insight into the hub structure of the causal human connectome, but it is not without limitations. One potential limitation lies in the use of a relatively coarse (n = 12) RSN partition for summarizing cortical hubs (Ji et al., 2019). However, the use of a published RSN partition facilitates interpretation of results, as the results of higher-dimensional (e.g., ICA-based) RSN partitions are often difficult to interpret and require abstraction via multidimensional statistics to summarize. In the Supplement we report the consensus network structure of these same data using a larger number of RSNs [the 22 neuroanatomical regions reported in the supplement of (Glasser et al., 2016a)]. The consensus structure of connectivity between these 22 regions also shows an orderly progression of information from visual sensory regions to the dorsal and ventral visual streams, through the parietal association cortex, and into the motor and prefrontal cortex. The organization of this more granular causal network also aligns well with recent perspectives on the hierarchical organization of the prefrontal cortex (Badre and Nee, 2018). An additional limitation of these results is that our method for calculating causal connectivity is unable to discover two-cycles (direct feedback cycles, where A->B and B->A). While two methods have been proposed for fMRI connectivity that are theoretically capable of recovering two-cycles (Sanchez-Romero et al., 2019), these methods have never been used in an applied research context, and often perform worse than the skew-based orientation method we use . Notably, the RSkew orientation method we adopt for the GANGO framework can discover 3-cycle or greater feedback loops, so only direct feedback loops remain unmeasured.
Nevertheless, as methods for more accurately assessing feedback cycles from fMRI are developed, the framework we propose in the current investigation could be expanded to include such methods.

Conclusion
Using a causal discovery machine learning framework, we demonstrate that the most centrally connected hubs of the cortical connectome are situated in the frontoparietal, dorsal attention and cinguloopercular networks. In particular, the causal human connectome highlights high connectivity of the frontoparietal network with all other higher cognitive RSNs. The discovered hub role of the frontoparietal network in the causal human connectome is especially attractive, as brain-based therapies for psychiatric conditions typically impact or directly stimulate nodes in the frontoparietal network (Belsher et al., 2021;Ferrarelli and Phillips, 2021;Fitzgerald, 2021;Song et al., 2019;Voigt et al., 2021;Zhang et al., 2021;Zilverstand et al., 2017Zilverstand et al., , 2016. Previously, we even demonstrated that connectivity in the frontoparietal network has downstream causal effects on the severity of alcohol use disorder (Rawls et al., 2021). As it is applied on a single-subject basis, the GANGO method could potentially even enable individualized causal connectivity-based neuromodulation targeting. Thus, the current study sets the stage for future applications of data-driven causal connectivity applications in psychiatry.