Brain effective connectome based on fMRI and DTI data: Bayesian causal learning and assessment

Neuroscientific studies aim to find an accurate and reliable brain Effective Connectome (EC). Although current EC discovery methods have contributed to our understanding of brain organization, their performances are severely constrained by the short sample size and poor temporal resolution of fMRI data, and high dimensionality of the brain connectome. By leveraging the DTI data as prior knowledge, we introduce two Bayesian causal discovery frameworks -the Bayesian GOLEM (BGOLEM) and Bayesian FGES (BFGES) methods- that offer significantly more accurate and reliable ECs and address the shortcomings of the existing causal discovery methods in discovering ECs based on only fMRI data. Moreover, to numerically assess the improvement in the accuracy of ECs with our method on empirical data, we introduce the Pseudo False Discovery Rate (PFDR) as a new computational accuracy metric for causal discovery in the brain. Through a series of simulation studies on synthetic and hybrid data (combining DTI from the Human Connectome Project (HCP) subjects and synthetic fMRI), we demonstrate the effectiveness of our proposed methods and the reliability of the introduced metric in discovering ECs. By employing the PFDR metric, we show that our Bayesian methods lead to significantly more accurate results compared to the traditional methods when applied to the Human Connectome Project (HCP) data. Additionally, we measure the reproducibility of discovered ECs using the Rogers-Tanimoto index for test-retest data and show that our Bayesian methods provide significantly more reliable ECs than traditional methods. Overall, our study’s numerical and visual results highlight the potential for these frameworks to significantly advance our understanding of brain functionality.


Introduction
Causal connectivity, also known as the Effective Connectome (EC), refers to the causal interaction between different regions within the brain [1].The EC focuses on understanding how neural activity in one brain region influences the activity in another region, thus revealing the flow of information and causal interactions between different parts of the brain.Causal discovery of brains mostly focuses on deriving the network techniques like stimulating the specific regions or build models based on physiological understandings, such as Dynamic Causal Modeling (DCM) [2].However, these methods are limited in their ability to perturb only a small number of dimensions among the multitude of regions involved, rendering this type of analysis impractical for high-dimensional ones.As a result, studies rely on functional connectivity, which is derived based on the correlation of neuronal activity of the brain regions, such as [3][4][5][6].However, the assumption that correlation implies causation is an erroneous approach and the main challenge lies in the disparity between the demand for high-dimensional causal statements and the limitations of providing low-dimensional causal statements.
The casual discovery of brain networks faces two main challenges: one with the fundamental impotence of methods and the other with the quality of data.While various computational methods, such as Bayesian network modeling [7], Granger Causality (GC) [8], and Structural Equation Models (SEM) [9] have been proposed, they tend to fall short from the perspective of accuracy and reliability.In particular, the Graphical approaches, or the Directed Acyclic Graph (DAG) frameworks, are shown to be more accurate in identifying the underlying causal mechanism based on observational and experimental data [10].As the graphical approaches, the Fast Greedy Equivalent Search (FGES) [11], NOTEARS [12], and Gradient-based Optimization of DAG-penalized Likelihood for learning linEar DAG Models (GOLEM) [13] have higher accuracy in large-scale networks, which makes them more appropriate in discovering EC.Assuming a static structure of EC, previous studies such as [14] and [15] demonstrate the effectiveness of the FGES and NOTEARS methods in discovering ECs.However, extracting causality in high-dimensional graphs and a limited number of samples still remains the main computational challenge.As shown in [13], as the network size increases and the sample size becomes limited, the accuracy of the derived graph significantly decreases.Aside from computational challenges, methods' performance can be affected by practical challenges like the low temporal resolution of fMRI, leading to results that are contaminated with fundamental ambiguities such as low reproducibility in discovered connectomes and less accurate results [16,17].In [18], it is shown that prior information that favors sparsity (such as DTI signals [19]) can significantly enhance both the accuracy and reproducibility of discovered graphs when a limited amount of data (such as fMRI signals) is available and the dimension of networks grows (such as brain network).Therefore, developing Bayesian causal methods based on multimodal data may help overcome the abovementioned challenges.
The DTI technique tracks nerve fibers as they connect potentially functionally associated brain regions.The recent research on multimodal data indicates that these nerve fibers (structural connectivity) constrain the strength and persistence of effective/functional connectivities [20][21][22][23][24][25][26].However, the results in [27,28] suggest that while there is some degree of overlap between functional and structural connectivity, there are also cases in which they do not align perfectly,i.e., there can be edges in functional connectivity matrix that are not in the binarized structural connectome (SC B ).The main reason for this misalignment is that functional connectivity refers to the statistical dependencies between different brain regions, suggesting that the structural connectome does not serve as reliable prior knowledge when exploring functional connectivity.In contrast, in the context of causal connectivity, there is an alignment between causal connectivity and structural connectivity, that is, the edges in the EC represent a subset of the edges that are present in the SC B [29].It is worth noting that a causal connection between two regions cannot be discovered without the mediation of another variable, assuming that the corresponding edge in SC B does not exist [30].Additionally, the subcortical regions can potentially act as mediating or confounding variables of cortical regions.As a result, including the subcortical regions in the causal discovery process is essential to ensure the satisfaction of the causal March 12, 2024 2/23 sufficiency assumption which is presented in [10].Studies such as [31] have shown that developing Bayesian frameworks for discovering ECs and leveraging multimodal data, such as DTI, as prior knowledge improves inference on effective connectivity, compared with unimodal studies.This research path has garnered significant attention in neuroscience and shows promise in understanding brain organization.Notably, as a causal modeling approach, the DCM method benefits significantly from extracting effective connections with priors on structural connectivity.In [32] and [33], the ECs for patients with depression and autism, respectively, are extracted by focusing on the structurally connected regions.Specifically, they apply a linearized version of the DCM to only the structurally connected regions, SC B .In [34,35], the authors presented a Bayesian DCM method using probabilistic structural connectivity, SC P , as prior information.The derived results in these studies are compared with fMRI-based models to illustrate the improvement of the accuracy of the EC model when SC P is employed as prior information.Moreover, causality in the brain is investigated in [36] with the Bayesian partial correlation method, in which the DTI data is approximated with the G-Wishart distribution as the conjugate prior.Likewise, in [37], the authors investigate causality by minimizing the reconstruction error of an autoregressive model constrained by the structural connectivity prior.
In this paper, we introduce two new Bayesian causal frameworks, i.e., Bayesian GOLEM (BGOLEM) and Bayesian FGES (BFGES), using DTI data as prior knowledge of the EC discovery to address the mentioned challenges of discovering EC when relying solely on fMRI data.Then, we show that these frameworks provide more accurate and reliable (reproducible) ECs that highlight the potential for these frameworks to advance our understanding of brain function and organization significantly.To compare the accuracy of different methods, we present the Pseudo False Discovery Rate (PFDR) metric based on the False Discovery Rate (FDR) as a computational metric in evaluating the accuracy of discovered ECs.Then, we measure the reliability of discovered ECs using the Rogers-Tanimoto index for test-retest data.The main contributions of our work are thus three-folded: • First, we introduce the BGOLEM and BFGES methods as Bayesian causal frameworks to improve the accuracy and reliability of discovering brain effective connectomes while addressing mentioned practical and computational challenges.These methods are Bayesian extensions of the well-established GOLEM and FGES methods.
• Second, we introduce the concept of the Pseudo FDR (PFDR) metric to illustrate the effectiveness of Bayesian causal frameworks, adjust each method's hyperparameters, and compare the accuracy of different methods.Our findings pinpoint the trustworthiness of PFDR as a measurable counterpart of FDR on real data.
• Third, we demonstrate that our Bayesian methods for discovering ECs achieve higher accuracy, as computed using the PFDR metric, and greater reliability, as assessed by the reproducibility metric known as the Roger-Tanimoto index.
The notable PFDR values computed for ECs discovered with non-Bayesian methods with unimodal data depict the weakness of these methods in discovering EC and emphasize the importance of employing Bayesian frameworks with multimodal data that lead to significantly accurate ECs.
The rest of the paper is organized as follows.The Materials and methods section presents our Bayesian frameworks, BGOLEM and BFGES methods, and our proposed computational accuracy measure, the PFDR.In the Results section, we utilize the methods by comparing different types of errors.Then, we employ Bayesian causal frameworks with HCP data and assess the accuracy and reliability of the discovered ECs through various means.Finally, we discuss the results and limitations of our methods.

Materials and methods
In this section, we begin by describing how to compute the probabilistic structural connectome.Then, we introduce our Bayesian causal frameworks, BGOLEM and BFGES, which are developed based on the GOLEM and FGES methods, respectively.In the BGOLEM method, we incorporate prior knowledge into the optimization process of the GOLEM method by masking it on the optimization score.The GOLEM method, which is based on the Gradient-based Optimization of DAG-penalized Likelihood for learning linear DAG Models, provides a foundation for our Bayesian framework.Similarly, in the BFGES method, we integrate prior knowledge with the Fast Greedy Equivalent Search (FGES) algorithm.In our Bayesian framework, we introduce a modified Bayesian information criterion (BIC) that takes into account the prior knowledge.Then, due to the absence of a gold standard for assessing the accuracy of effective connectivity (EC) methods, we introduce the concept of the PFDR metric as a means to evaluate the performance of the BGOLEM, BFGES, GOLEM, and FGES methods.The PFDR metric is conceptualized with an eye on the idea from the DTI technique and the FDR metric.In the data subsection, we begin by discussing the Causal sufficiency assumption and its relation with our choice for parcellation, explaining its significance and implications for our study.We then outline the reasons for selecting a specific parcellation approach, highlighting its relevance and benefits for our research.Following that, we present the HCP data used in our study, providing details about the data acquisition and preprocessing procedures.Next, we describe the process of generating synthetic and hybrid data, elaborating on the methodologies and considerations involved.We explain how these data sets complement the empirical data

Probabilistic structural connectome
The prior knowledge in the developed Bayesian methods is computed from DTI data.As a result, in order to compute the probabilistic SC, SC P , we use the following formula where ρ ij represents the streamlines from seed i to the target region j that is normalized according to the area of the j th region [34] and n is the number of regions in the brain atlas.

Bayesian causal discovery frameworks
Consider a parameterized Bayesian network model with n nodes (B), defined by (G, θ), where, G = (V, E) is a DAG, V is a set of nodes (U = {X 1 , ..., X n }), E is a set of directed edges (causal relations), and θ is a set of parameters that specify all conditional probability distributions.A parameterized Bayesian network B represents and factors a joint distribution over U according to the structure where P a G i is the set of parents of X i in G, and θ i ⊂ θ is the subset of parameters that define the conditional probability of node X i ∈ V given its parents in G.

Proposed Bayesian GOLEM method
The GOLEM algorithm finds a linear DAG model that equivalently represents a set of linear structural equations.In GOLEM, P B in equation 2 is generated by the linear DAG, X = W X + N , where, X = [X 1 , ..., X n ] T , W is a n × n weighted adjacency matrix, and N = [N 1 , ..., N n ] T is an independent noise vector.In [12], W matrix is found by optimizing a score function F B (W, X), subject to the structure G, for a given samples X where G(W ) is the n-node graph induced by the weighted adjacency matrix W and F B (.) : R n×n → R is the score function.l(W, X) is the maximum likelihood, R sparse (W ) is a penalty term that favors sparsity, i.e., having fewer edges.The hard optimization problem in equation 3 is relaxed into a soft optimization problem presented in [13] as follows min March 12, 2024 5/23 where R DAG (W ) is a penalty term that favors DAGness of W .The penalty term for encouraging sparsity is defined as λ GOLEM ∥W ∥.The details of the algorithm and proofs are found in [13].
BGOLEM According to [13], the GOLEM method results in lower Structural Hamming Distance (SHD) and FDR values for the sparse structures.However, for a limited sample size, the accuracy of GOLEM significantly decreases with a higher number of edges and nodes.To cope with this deficiency, we define the BGOLEM algorithm by suggesting the following sparsity penalty term in 4, to include the prior probability where Q(.) : R n×n × R n×n → R n×n is a smooth element-wise monotonic function.This function acts element-wise on the input matrices, such that if a specific element in the P (G) matrix moves in one direction, the corresponding element of the W moves in the same direction.Such a sparsity penalty term inhibits/excites an element of the W matrix that has a lower/higher prior probability.In this paper, we have considered where ⊙ is the element-wise matrix product.The choice of the Q(.) function can change the incorporation of prior information in 5. Employing the log(.)function equalizes the skewed distributions, in a highly right-skewed distribution such as the SC P , where most of the values are close to zero, the logarithmic transformation can be beneficial in equalizing the distribution.The impact of extreme values is reduced, which can help bring the distribution closer to a normal or symmetric shape.In the context of incorporating prior knowledge into the final graph, the logarithmic transformation results in a more conservative approach.This conservative incorporation of prior knowledge results in a more balanced representation in the final graph.Moreover, the logarithmic transformation can help reduce the impact of noise on the data.
Compressing extreme values mitigates the influence of outliers and makes patterns in the data more apparent.To find a proper λ BGOLEM value, an exhaustive search can be performed to find the optimum λ BGOLEM value with the lowest false discovery rate.However, for fMRI data, since there is no access to the ground truth, we use the PFDR metric that is presented in The PFDR metric section to find the optimum hyperparameters, i.e.,λ BGOLEM , value.

Proposed Bayesian FGES method
The FGES causal discovery method presented in [11], has two main parts that are: how the algorithm operates and how to define the scoring criteria.This algorithm has two steps and starts with an empty graph.In the first step, the forward phase, scores of all alternative one-edge additions to the graph are computed.The edge that corresponds with the best score is added to the graph.This process continues until no more improvement is achieved by adding a single edge.In the second step, the backward stage, the graph is pruned backward elimination through step-wise single-edge deletion.This procedure continues until any single-edge deletion results in a decline in the score.
The FGES uses a decomposable score where there is no need to compute the score of the entire graph in each iteration.Only the scores of the nodes with changing parents change.This is the main difference between the FGES and GES methods [11,38].The March 12, 2024 6/23 decomposable scoring criteria S B [11], is defined by where D is a set of observed data and s(X i , P a G i ) is the score corresponding to a node X i and its parents P a G i .The Bayesian score function S B using the relative log posterior of G is where p(G) is the prior probability of G, and p(D|G) is the conditional likelihood of D [11] In particular, in [39], the authors show that equation 8 for curved exponential models can be approximated using Laplace's method for integrals, yielding where θ denotes the maximum-likelihood estimate for the network parameters, d denotes the dimension (i.e., number of free parameters) of G, and m is the number of records in D. In [11], the decomposable score for the FGES method is defined by the sum of the second and third terms in equation 9 that is also known as the BIC, i.e., where λ FGES is the hyperparameter of the BIC score.
BFGES Similar to the GOLEM method, the FGES method has difficulties discovering DAGs with a higher number of nodes and edges, and a limited number of samples, as illustrated in [12].While the log p(G) term is neglected in the BIC scoring 10, this term can play a dominant role in model selection of large-scale graphs [18].In our proposed Bayesian FGES algorithm, considering the equations 7 and 8, we derive the decomposable version of equation 9 as follows where d i denotes the number of parameters for the structure of X i and its parents.The procedure of deriving equation 11 is quite similar to that of equation 10, but without ignoring the prior probability.If we assume a uniform prior for the existence of an edge, all three probabilities of X i ← X j , X i → X j and X i . . .X j (X i and X j are not connected) are equal to 1 3 .Since the probability of X i not being a parent of "X j is equal to 2  3 ", then p(G) = 1 3 −N G , where N G an n are the number of edges and number of nodes in graph G, respectively.Similar to the FGES algorithm, in the BFGES algorithm, we start with an empty graph, which means that p(G) is the probability of a graph G with no edges.As a result, the algorithm starts with March 12, 2024 7/23 where p(X i ↔ X j ) = p(X i → X j ) + p(X i ← X j ) is the edge existence probability between nodes i and j, i ̸ = j.Consequently, effect of adding an edge from node i to node j is to subtract log 1 − p(X i ↔ X j ) from the score and add log 1 2 p(X i ↔ X j ) to the score where in this equation, S B (G ′ , D) is the score of graph G ′ that is derived from adding In the backward stage of the BFGES algorithm, the score of the derived graph is subtracted with log 1  2 p(X i ↔ X j ) and log 1 − p(X i → X j ) is added when the edge from node j to i is eliminated where In sum, the BFGES algorithm runs similarly to the FGES algorithm with the modified scoring 11 to incorporate prior knowledge in causal discovery.Similar to the BGOLEM method, for fMRI data, we use the PFDR metric (The PFDR metric section) to find the optimum λ BFGES value.

The PFDR metric
To measure the reliability of methods in discovering causality, False Discovery Rate (FDR) is a powerful full metric to assess the accuracy of discovered edges, which is defined as follows where F P and N oDE are the number of False Positives and the Total Number of Discovered Edges, in the true underlying EC, respectively [40].In the assessment of discovered ECs, there is no access to the true causal underlying graph.As a result, computing F P and consequently, FDR value is not possible.
In the DTI technique, the relative diffusivity of water in a voxel into directional components is quantified.The longest axis of the diffusion ellipsoid, which is estimated based on diffusion tensors, is used to track nerve fibers as they travel between potentially functionally associated brain regions [41].The absence of an edge in SC B indicates that there is no physical connection between two regions which implies that no effective connection is possible.On the other hand, the existence of an edge in this matrix implies a physical path between two regions, which opens the door for effective connectivity.Therefore, from the absence of an edge in this matrix, one can infer that the corresponding element of the EC must be zero.The presence of an edge in SC B allows for corresponding non-zero elements in EC.We exploit DTI-based SC B information and propose the Pseudo FDR (PFDR) metric, which is defined as follows where computation, as a result, F P ′ ≤ F P .Note that PFDR and FDR can be equivalent in the case the undirected causal graph of the underlying brain's causal mechanism exactly matches the structural graph.Contrary to the FDR metric, computing the PFDR metric is practicable, while we do not know the true underlying graph.This metric enables us to measure and compare the accuracy of the ECs discovered with causal discovery methods.The PFDR can be used to adjust the hyperparameters of causal discovery methods, as well.In the Results for the empirical data section, we experimentally explore the relationship between FDR and PFDR by deriving both values for hybrid data and computing their correlation.This analysis aims to emphasize that these two metrics can be used interchangeably.Details and discussions on Pseudo FDR can be found in Appendix 1.

Data
Casual sufficiency and selected Parcellation The causal sufficiency assumption is the concept that addresses the question that whether any extra variable need to be considered in the process of the causal discovery process of a causal model [10].This concept refers to the condition where, if a set of variables X are causally sufficient, it should not omit any hidden common cause C / ∈ X that has the ability to influence more than one variable within X [42].Since the subcortical regions receive information flow originating from the brain stem, they have the potential to act as confounding variables to the cortical regions and may play a role as confounding variables in their causal relationships.As a result, to ensure that the causal sufficiency assumption is satisfied, the subcortical regions are included in our analysis.Furthermore, consider the scenario where region A comprises sectors a 1 and a 2 , while region B consists of sectors b 1 and b 2 .In this case, sector a 1 might exert a causal influence on sector b 1 , while sector b 2 could have a causal influence on sector a 2 .Such interdependencies between sectors can pose challenges during the causal discovery process and potentially disrupt the identification of causal relationships.To mitigate these challenges, it becomes crucial to select an atlas that encompasses a rich array of brain regions.By utilizing an atlas with a greater number of distinct regions, we increase the likelihood of capturing the individual sectors involved in the causal relationships accurately.Correspondingly, we employ the Destrieux atlas [43] that satisfies our concerns.The List of ROIs of the Destrieux atlas is presented in Appendix 2.

Data acquisition
In this study, we used DWI and resting-state fMRI data from unrelated subjects of the "HCP1200" data set (March 2017 data release of healthy adults aged [22][23][24][25][26][27][28][29][30][31][32][33][34][35] [44].The HCP data sets were acquired using protocols approved by the Washington University institutional review board, and written informed consent was obtained from all subjects.All the authors have been approved by HCP to use this data. DTI data preprocessing DWI images were acquired using a 3T 'Connectome Skyra', provided with a Siemens SC72 gradient coil and stronger gradient power supply with maximum gradient amplitude (Gmax) of 100 mT /m (initially 70 mT /m and 84 mT /m in the pilot phase) sequence and 90 gradient directions equally distributed over 3 shells (b-values 1000, 2000, 3000 mm/s 2 ) with 1.25mm isotropic voxels [45].The preprocessing steps included distortion and motion correction, brain extraction, and eddy correction for diffusion images, as well as registration of structural and diffusion data to each other in native space [46].In this study, we employed probabilistic tractography that considers multiple possible pathways at each voxel to find any possible structural connection.These preprocessing steps were already performed in the downloaded package from HCP web page.The MRtrix3 software is used to perform March 12, 2024 9/23 whole-brain probabilistic tractography, which involved using the 5TTgen function with the FSL option to segment subcortical nuclei in T1-weighted images.The segmented subcortical ROIs are then combined with cortical ROIs from the Desikan-Killiany atlas that are segmented by Freesurfer during the preprocessing steps in the HCP pipeline.The msmt 5 tt algorithm is used to estimate the response function for diffusion images, and the iFOD2 algorithm is used with the multi-tissue multi-shell -msmt c sd flag to generate FOD (fiber orientation density) images.The tckgen algorithm is then used to generate 50 million streamlines from the FOD images with a maximum tract length of 250 mm to obtain whole-brain probabilistic tractography.Finally, the -tck2connectome function is employed to find the number of streamlines that started from each ROI and ended in all other ROIs to determine the structural connectivity between all nodes.

Synthetic and hybrid data
The synthetic data, the first group of data, is generated using a random graph generated from the scale-free (SF) model, as detailed in [12].This process involves generating graphs with a degree of 7 and varying the number of nodes to be 50, 75, 100, 125, and 150, along with 300 time points.For each number of nodes, a total of 30 datasets are generated.This group of data is employed to compare the results and effectiveness of our Bayesian methods with the results in [13].
To generate hybrid data, the second group of data, two main steps are involved.The first step is to generate random causal graphs based on the DTI data from the HCP.These causal graphs represent the underlying DAG that defines the causal connectivity between different brain regions.Once the causal graphs are established, the second step is to generate fMRI data based on these causal graphs and the hemodynamic response function.This involves simulating the neural activity within the brain regions defined by each causal graph and then convolving this neural activity with the hemodynamic response function to generate the corresponding fMRI signals.We create 50 causal graphs with 164 nodes and the final hybrid data consists of 1200 time points.The hybrid data is used to show the validity of the PFDR metric.
More in-depth information on the generation of Hybrid data can be found in Appendix 1.

Empirical data
The HCP MRI data acquisition protocols and procedures have previously been described in full detail [47].We used the minimally preprocessed images of resting state fMRI, 'rsfMRI', and Diffusion Tensor Imaging (DTI), 'dMRI', for 50 unrelated subjects that were provided by the HCP S1200 Release.fMRI resting-state runs (HCP file name: rfMRI-REST1 and rfMRI-REST2) were acquired in four separate sessions on two different days, with two different acquisitions (left to right, or LR, and right to left, or RL) per day [48].To eliminate any potential systematic bias or artifacts that may merge due to the acquisition process, one can average the two sessions of each day [49] or average the results of applying any connectome discovery methods to two sessions [50].In our study, for each subject, the mean of the first day's two sessions is used as test data (Mean 1), while the second day's two sessions are used as retest data (Mean 2).Averaging allows for a comprehensive analysis that incorporates data from both acquisition directions, enhancing the reliability and validity of the results.Note that in both empirical and hybrid data, we use DTI data of the same subjects.

Results
Applying the developed Bayesian methods, the probabilistic SC P is computed using equation 1.To obtain the binary structural connectome (SC B ) that is used in the PFDR metric, in the first step, the number of tracts for each subject is binarized using a threshold of 50% of stream counts.In the second step, majority voting is applied to the results of the first step for 50 subjects to compute the final binary SC, SC B .It is worth noting that studies such as [19] have shown that it is reasonable to exclude at least 80% of the edges when deriving the binary SC B .However, in our study, we have opted to use a 50% threshold for binarization.This choice is motivated by the significance we attribute to zero values in the SC B used in the PFDR metric, as well as our intention to be more conservative in excluding edges that are not consistently connected across our subjects.By employing majority voting on zero values, we aim to ensure that the derived SC B is as conservative as possible in capturing the absence of connections.It is worth noting that extracted SC B has 33% zero elements.
In the following, to highlight the effectiveness of our Bayesian methods, we apply the GOLEM and FGES methods along with the BGOLEM and BFGES on the synthetic data (SF7) and show the results in terms of the FDR (False Discovery Rate) and the Percentage of Total Errors.The Percentage of Total Errors is defined as follows, False Positives + False Negatives Total Number of Discovered Edges (17) Then, to assess the reliability and validity of the proposed PFDR metric, we compute the correlation of the FDR and PFDR values of the ECs discovered by applying FGES, BFGES, GOLEM, and BGOLEM to the hybrid data.To apply the Bayesian causal frameworks to the empirical data, we utilize the SC P of each subject as prior knowledge for the EC discovery process of the same subject.Results for these data using all four methods are obtained by adjusting their hyperparameters based on the minimum PFDR values observed in the test data.In comparing the accuracy of the Bayesian method with non-Bayesian methods, to avoid any biases, the PFDR for each subject is computed using the SC B , which is obtained by excluding the same subject from a group of 50 subjects.In the end, to test the reliability of the ECs discovered using our Bayesian methods is examined using the Roger-Tanimoto index, and the results are compared with those obtained from the ECs discovered using the FGES and GOLEM methods.Visual comparisons of the ECs discovered with different methods are also provided to facilitate a comprehensive analysis.

Results for the synthetic and hybrid data
In this subsection, we present a comparative analysis of the results and effectiveness of our Bayesian methods in contrast to the outcomes obtained using the traditional method with synthetic data.For Hybrid data, both the FDR and PFDR values are computable.Subsequently, by comparing the PFDR and FDR values of applying both Bayesian and non-Bayesian methods to hybrid data, we evaluate the trustworthiness and validity of the PFDR metric.As a result, the optimal values of λ FGES , λ GOLEM , λ BFGES and λ BGOLEM are determined based on minimizing the Percentage of Total Errors, as ground truths are available for these two groups of data.
The results obtained from applying the FGES and BFGES methods, along with the GOLEM and BGOLEM methods, on synthetic data with varying numbers of nodes are depicted in Fig 2 and 3.

Results for the empirical data
In this section, we present and compare the results of applying our methods and traditional methods to empirical data, both numerically and visually.In Table 1, we compare the PFDR values for the Mean 1 (test) and Mean 2 (retest) data when 200 edges are selected.To assess the significance of employing our Bayesian methods compared to the non-Bayesian method in terms of accuracy, we conducted the Paired T-test on the computed PFDR values.This test is selected to be performed on the PFDR values due to appropriation of this test for paired observations that the differences between pairs are normally distributed.The p-value obtained for comparing the PFDR values of BFGES and FGES is less than 1e-3, indicating a significant difference between these two methods.Similarly, the p-value for comparing the PFDR values of BGOLEM and GOLEM is also less than 1e-3, indicating a significant difference between these two methods as well.The PFDR metric compares the accuracy of each method in parts of the discovered ECs that are zero in the SC B and there is still no ground truth about the parts of ECs that are one in the SC B ., i.e., the PFDR metric is able to identify the part of the false positive in the ECs.We used the proportion test to compare the percentage of significantly different edges in ECs.In doing so, we searched for edges that have significantly different results by using different methods, namely FGES and GOLEM, in the population of 50 subjects.Our results show that 40.3% of EC edges that are zero in SC B are significantly different between these two methods, whereas this number is reduced to 11.2% for the rest of the edges, suggesting that edges that are zero in SC B are more sensitive to the used method.
To further investigate the effectiveness of the Bayesian causal frameworks on the reliability of discovered directed (asymmetric) and undirected (symmetric) ECs, we use the Rogers-Tanimoto index [51] that evaluates the reproducibility of each edge within the derived ECs.Note that, in undirected ECs, the connectome is primarily determined by the presence or absence of an edge, leading to a symmetric connectome.However, in directed ECs, both the existence or absence of edges and the direction of edges are discovered, leading to an asymmetric connectome.This index calculates the dissimilarity between the edges of the ECs obtained from the test and retest data of 50 subjects.By calculating this index for all the edges, we gain insights into the degree of alignment between the extracted ECs in the test and retest data.The resulting index value for each edge ranges from 0 to 1, where 0 indicates a perfect agreement and no dissimilarity of an edge in two sets of derived ECs, while 1 denotes no agreement and complete dissimilarity between an edge.In  To assess the significance of employing our Bayesian methods compared to the non-Bayesian method in terms of reliability, we conducted the Wilcoxon signed-rank test on the computed Roger-Tanimoto values.This test is performed due to its appropriation for the comparison of two sets of data measured on an ordinal scale that are highly skewed.The p-value of this test obtained for comparing the Roger-Tanimoto values of BFGES and FGES is extremely small, with a value lower than 1e-3.Similarly, the p-value for comparing the Roger-Tanimoto values of BGOLEM and GOLEM is less than 1e-3.These results indicate that the observed differences in the Roger-Tanimoto values between BFGES and FGES, as well as between BGOLEM and GOLEM, are highly significant.The Rogers-Tanimoto index for undirected ECs is illustrated and discussed

Discussion
In this paper, first, the Bayesian causal frameworks are introduced in the Proposed Bayesian GOLEM method and Proposed Bayesian FGES method subsections.Then, the PFDR metric is presented in The PFDR metric subsection as a computable alternative for the FDR metric to assess the accuracy of causal discovery methods.Fig 2 and 3 show that the mean FDR and the mean Percent of Total Error are improved with the employment of Bayesian causal frameworks for SF7 data with different number of nodes.The results of applying the FGES and GOLEM methods to this data are consistent with the results in [13].Fig 4, shows the relationship between FDR and PFDR using hybrid fMRI data with 164 nodes.The correlation coefficient between the PFDR and FDR of ECs obtained using the BGOLEM, GOLEM, BFGES, and FGES methods are close to one, i.e., the PFDR and FDR values are highly correlated in all of the methods.These results depict that the PFDR metric can be used as an accuracy metric instead of the FDR metric for discovering ECs where ground truth data is not available.It is worth noting that in the traditional method for evaluating the accuracy of causal discovery methods in extracting ECs, the first step involves collecting the differences in ECs discovered by applying a specific causal discovery method to both patient and normal subjects' data.Then, the aim is to assess the consistency between the discovered differences and the findings of neuroscientific studies.This type of evaluation is used in studies such as [15,32,33].While this approach can be informative about some causal interactions in brain networks, the agreement/disagreement of small parts of ECs of patients and normal subjects does not necessarily imply the accuracy/inaccuracy of a causal method.For instance, an accurate functional connectome discovery method can identify a difference between two groups of data (patient and normal), however, identifying this difference does not imply that the derived connectomes are presenting the causal network of the brain.
Fig 5 illustrates the structural and functional connectomes along with ECs discovered with Bayesian and non-Bayesian methods in various aspects.This figure investigates the visual symmetry of connections in the two hemispheres and the edges between them.Comparison between the functional connectome derived from correlation values and the structural connectomes reveals that several regions are functionally highly correlated despite lacking structural connections.This functional connectome resembles those in [52].Moreover, studies, such as [53][54][55], have extensively investigated brain asymmetry and provided insights into the functional differences between the two hemispheres.However, the functional connectome depicted in Fig 5 exhibits high symmetry between the right and left hemispheres, which emphasizes of shortcomings of this method in understanding brain mechanisms.As a result, this contrast between the functional and structural connectomes underlines the opposition between the two, suggesting that the structural connectome can not serve as reliable prior knowledge when exploring functional connectivity.This is expected due to the existence of confounding nodes and the fact that functional connectivity refers to the statistical dependencies between different brain regions.These findings highlight the results presented in [27] and [28].
Fig 5 helps to visually evaluate the potential of our results.From a structural perspective, it is well established that the two hemispheres of the brain are connected through commissural nerve tracts such as the Corpus callosum, which bridges the left and right hemispheres to share information.According to the structural connectome presented in Fig 5, with the exclusion of subcortical regions of the brain, there is no structural direct edge between the two hemispheres.As a consequence, although we observe a high correlation between regions of two brain hemispheres, due to the satisfaction of the causal sufficiency assumption, it is not possible to establish a direct causal relationship between hemispheres.In this figure, the ECs discovered through the FGES and GOLEM methods show connected nodes between the two hemispheres, with many edges having a low probability of existence in each hemisphere.Moreover, while the ECs discovered through effective methods can not be symmetrical between the two hemispheres, the ECs discovered with traditional causal discovery methods are highly symmetric, as it is shown in Fig 5.This emphasizes the limitations of discovering EC using state-of-the-art causal discovery methods like FGES and GOLEM, particularly when relying solely on fMRI data.These challenges can arise due to computational limitations (the weak performance of state-of-the-art methods in dealing with high-dimensional networks with short sample sizes) as well as practical constraints ( poor temporal resolution of fMRI data).In contrast, the ECs discovered with Bayesian frameworks have a low number of edges between the two hemispheres, with a greater proportion of edges having a higher probability of existence in each hemisphere.These visual observations, along with the requirement that EC edges should be a subset of the edges in the SC, highlight the significance of employing Bayesian approaches in the investigation of EC.
To numerically evaluate the effectiveness of our method in the accuracy of derived ECs, we employed the PFDR metric and the results are presented in Fig 6 .Utilizing the GOLEM and FGES methods with empirical data shows that the GOLEM method has a lower mean PFDR value than that of the FGES method, which implies that the GOLEM method is more accurate in discovering EC (Fig 6 [14] and [15].Nevertheless, there is still no measure to assess the accuracy of causal discovery methods in determining the presence or absence of an edge in the EC, that is, one in the SC B .The extracted ECs with different methods vary across subjects; however, this variance is more pronounced in areas of the EC whose corresponding edges are zero in SC B .To support this statement, we used the proportion test to compare the percentage of significantly different edges in ECs when 200 edges are selected.In doing so, we searched for edges that had significantly different results by using different methods, namely FGES and GOLEM, in the population of 50 subjects.Our results show that 40.3% of edges that are not connected in SC B are significantly different between these two methods, whereas this number is reduced to 11.2% for the rest of the edges, suggesting that edges that are not connected in SC B are more sensitive to the used method.Although we did this analysis for 200 edges in the ECs, similar results were obtained for other numbers of edges. The reliability of EC discovery methods is numerically assessed using the Roger-Tanimoto index as a measure of reproducibility.Fig 7 demonstrates that both directed and undirected edges between the two hemispheres exhibit greater similarity in the Bayesian methods, as indicated by their Roger-Tanimoto values.Furthermore, Fig 8 illustrates a decrease in the Roger-Tanimoto index for the ECs when employing the Bayesian versions of the methods.The Bayesian causal frameworks exhibit significantly higher reproducibility (lower dissimilarities) in discovering directed ECs compared to the non-Bayesian FGES and GOLEM methods.Moreover, the BGOLEM method exhibits higher reproducibility compared to the BFGES method.The Roger-Tanimoto index values for the undirected ECs decrease when the Bayesian causal frameworks are employed, as shown in the appendix.According to the p-values of the Wilcoxon signed-rank tests, utilization of our Bayesian methods significantly improves the reliability of the discovered ECs compared to those derived in [14] and [15] with the most state of the art methods, the FGES and GOLEM.It is important to note that the observed differences in the derived ECs using our Bayesian frameworks could be attributed to the intrinsic variations in brain connectivities among different subjects and not necessarily the weakness of causal discovery methods, as demonstrated in [56].

Limitations and promising aspects of future research
Future research needs to delve more deeply into four key limiting and challenging aspects.First, according to the findings of this paper, some errors that have already been demonstrated in the presented methods are caused by the limitations of fMRI in presenting the brain activities and causal interactions in the brain [57,58].As a result, in addition to DTI data, combining fMRI data with other neuroimaging techniques with a higher temporal resolution (such as EEG or MEG) can be beneficial in enhancing the accuracy of brain networks, as it is shown in [59].In this paper, we have shown the improvement of the discovered ECs by employing the Bayesian causal frameworks on both hybrid and real-world data.The degree of incorporation of the structural data in each method is optional, and various functions can be used to mask the prior information in each technique to discover a more accurate EC.However, there is no exact procedure for finding the optimal value that balances the impact of DTI and fMRI data on discovering ECs.This balancing value for fMRI and DTI fusion plays an essential role in the discovered ECs, and finding this value can be a second aspect of future research.
Moreover, considering multivariate Gaussian density distribution to estimate fMRI signals or assuming linear causal interaction between brain regions are fairly restrictive assumptions, as assumed in the FGES and GOLEM methods, respectively.As the third aspect of future works, since the activity in one neuronal population's gates connection strengths among others, and this procedure causes highly nonlinear information transformation between brain regions, model-free methods need to be developed to measure these interactions and discover the EC.
Most recent research investigates functional connectivity while considering time as one of the important parameters of their research, such as [60,61], which can provide information about a possibly low-dimensional and intrinsic manifold of brain data, as it is mentioned in [62].It is presented in [63] that discovering dynamic brain networks can lead to a revolution in our understanding of brain functionality and the future promises of this research path are demonstrated in [64].However, the existing research that is concerned with extracting effective connectivity from observational and interventional data is concerned with deriving static networks for brain EC and trying to interpret causal relations in the brain with them, and few research such as [65] have considered time as one the main component of EC discovery.As the fourth and most important aspect of future studies, the next generation of EC discovery from observational and interventional data must derive dynamic effective networks (dynamic effective connectomes) to exploit temporal and spatial information.We believe this research path can lead to an even more profound understanding of brain mechanisms than dynamic brain networks.

Conclusion
The ultimate goal of this paper is to develop more accurate and reliable methods for causal discovery in brain connectomes and to demonstrate their effectiveness both numerically and visually.We introduce two Bayesian causal discovery methods that leverage DTI as prior knowledge.We then demonstrate the effectiveness of our methods through a series of simulations on synthetic and hybrid data.To assess the accuracy of the derived effective connectomes with our methods, we introduce the PFDR metric as a computational accuracy assessment.We show that our methods produce more accurate effective connectomes.Furthermore, we use the Roger Tanimoto index as a reproducibility metric to demonstrate that the effective connectomes of test-retest data derived from our methods are more reliable than those derived from traditional methods.Our study emphasizes the potential of these frameworks to significantly advance our understanding of brain function and organization.

Fig 1 .
Fig 1. Steps in Bayesian causal learning and Assessment of effective connectomes.Deriving effective connectomes for both hybrid and empirical data with introduced Bayesian causal frameworks, computing the PFDR value and performing reliability tests.

Fig 2 Fig 2 .Fig 3 .
Fig 2. The mean FDR and the mean Percent of Total Errors for graphs with different number edges.(a): The mean FDR for the FGES and BFGES methods (b): The mean Percent of Total Errors for the FGES and BFGES methods

Fig 5 Fig 4 .Fig 5 .Fig 6 .
Fig 4. Dependency of the PFDR and FDR metrics for ECs derived with Bayesian and Non-Bayesian methods on hybrid data with 164 nodes.(a): The correlation coefficients (R) of PFDR and FDR values for the ECs of the GOLEM and BGOLEM are 98.9% and 99.9%, respectively (b): The correlation coefficients of PFDR and FDR values for FGES and BFGES are 80.9% and 97.7%, respectively.

Fig 7 ,Fig
we demonstrate the Rogers-Tanimoto values for all the edges of the ECs with 164 regions extracted using Bayesian and non-Bayesian methods for the test and retest data.The index is computed for both undirected and directed edges in this figure.Comparing the effect of employing Bayesian and non-Bayesian methods in the reproducibility of discovered asymmetric and symmetric ECs.

Fig 8
Fig 8 shows the numerical values of the Rogers-Tanimoto index for all the directed edges discovered with the GOLEM, BGOLEM, FGES, and BFGES for test and retest data.The median and interquartile range of the Rogers-Tanimoto index for directed ECs of the FGES method are 4% and [4, 11.5]%.These values for the directed ECs of the GOLEM methods are 4% and [4, 11.4]%.The median and interquartile range for the BFGES method are 4% and [0, 7.8]%, and for the BGOLEM method, these values are 0% and [0, 4]%.According to Fig 8, the computed values of this index for the ECs of the Bayesian and non-Bayesian methods are highly right-skewed.To assess the significance of employing our Bayesian methods compared to the non-Bayesian method in terms of reliability, we conducted the Wilcoxon signed-rank test on the computed Roger-Tanimoto values.This test is performed due to its appropriation for the comparison of two sets of data measured on an ordinal scale that are highly skewed.The p-value of this test obtained for comparing the Roger-Tanimoto values of BFGES and FGES is extremely small, with a value lower than 1e-3.Similarly, the p-value for comparing the Roger-Tanimoto values of BGOLEM and GOLEM is less than 1e-3.These results indicate that the observed differences in the Roger-Tanimoto values between BFGES and FGES, as well as between BGOLEM and GOLEM, are highly significant.The Rogers-Tanimoto index for undirected ECs is illustrated and discussed .a).Fig 6.b illustrates that both Bayesian methods decrease the PFDR values in EC discovery compared to the non-Bayesian methods.Considering the p-values of the Paired T-test, our numerical results imply that the proposed methods significantly decreased the errors in discovering accurate EC and enhanced the accuracy of discovering ECs in P ′ is the number of edges that are one in EC and zero in SC B and T P ′ is the number of edges that are one in both EC and SC B .From the definition of F P and F P ′ , false positives in PFDR computations are a subset of true false positives in FDR March 12, 2024 8/23

Table 1 .
Results of applying Bayesian methods on mean PFDR values for ECs with 200 edges for test and retest data The box plots of the Rogers-Tanimoto values of the FGES, BFGES, GOLEM, BGOLEM.The median and interquartile range of the Rogers-Tanimoto index for directed ECs of the FGES method are 4% and [4, 11.5]%.These values for the directed ECs of the GOLEM methods are 4% and[4, 11.4]%.The median and interquartile range for the BFGES method are 4% and [0, 7.8]%, and for the BGOLEM method, these values are 0% and [0, 4]%.