Integrative Bayesian Analysis of Brain Functional Networks Incorporating Anatomical Knowledge

Ixavier A. Higgins

Department of Biostatistics and Bioinformatics, Emory University, USA

Suprateek Kundu

Department of Biostatistics and Bioinformatics, Emory University, USA

Ying Guo^{1}

Department of Biostatistics and Bioinformatics, Emory University, USA

Abstract

Recently, there has been increased interest in fusing multimodal imaging to better understand brain organization by integrating information on both brain structure and function. In particular, incorporating anatomical knowledge leads to desirable outcomes such as increased accuracy in brain network estimates and greater reproducibility of topological features across scanning sessions. Despite the clear advantages, major challenges persist in integrative analyses including an incomplete understanding of the structure-function relationship and inaccuracies in mapping anatomical structures due to inherent deficiencies in existing imaging technology. This calls for the development of advanced network modeling tools that appropriately incorporate anatomical structure in constructing brain functional networks. We propose a hierarchical Bayesian Gaussian graphical modeling approach which models the brain functional networks via sparse precision matrices whose degree of edge-specific shrinkage is a random variable that is modeled using both anatomical structure and an independent baseline component. The proposed approach adaptively shrinks functional connections and flexibly identifies functional connections supported by structural connectivity knowledge. This enables robust brain network estimation even in the presence of mis-specified anatomical knowledge, while accommodating heterogeneity in the structure-function relationship. We implement the approach via an efficient optimization algorithm which yields maximum a posteriori estimates. Extensive numerical studies involving multiple functional network structures reveal the clear advantages of the proposed approach over competing methods in accurately estimating brain functional connectivity, even when the anatomical knowledge is mis-specified up to a certain degree. An application of the approach to data from the Philadelphia Neurodevelopmental Cohort (PNC) study reveals gender based connectivity differences across multiple age groups, and higher reproducibility in the estimation of network metrics compared to alternative methods.

Keywords: Adaptive shrinkage; brain networks; Gaussian graphical models; multimodal imaging; Philadelphia Neurodevelopmental Cohort; reproducibility.

## 1 Introduction

The human brain is an extremely complex organ responsible for all thought and bodily function. Various approaches have sought to explain the brain’s functionality as a result of neurotransmissions between individual neurons, reflected as the co-activation between voxels in brain images. Recently, there has been a rapid increase in research on brain connectome analysis focused on linking inter-regional dependencies to brain function. Methods for both structural connectivity (SC) and functional connectivity (FC) have seen increasing developments with the emergence of non-invasive technologies such as diffusion tensor imaging (DTI) and functional magnetic resonance imaging (fMRI). FC measures the temporal coherence in brain activity across two distinct brain regions, while SC approaches based on DTI reconstruct white matter pathways in the brain by measuring the diffusivity of water molecules in brain tissues. These two types of connectivity offer complementary and interdependent information about brain structure and function.

Despite strong evidence regarding the role of white matter fiber tracts in regulating FC (Damoiseaux and Greicius, 2009; Honey et al., 2010; Sporns, 2013) and considerable progress in separately estimating FC and SC, there have been comparatively limited advances in FC approaches which are guided by underlying anatomical knowledge. Incorporating anatomical knowledge in estimating FC is clearly desirable since it is expected to produce more accurate estimates of the network, which translates to greater reproducibility of the findings as illustrated via our fMRI data analysis. However, several considerations need to be taken into account, such as the complexity of the structure-function relationship (Hermundstad et al., 2013), heterogeneity in FC for a given SC strength, which presumably is attributed to the fact that FC is only partially dependent on SC (Damoiseaux and Greicius, 2009; Messé et al., 2014) and regulated by unobservable dynamics in underlying neuronal activity (Bressler and Tognoli, 2006).

Recently, Venkataraman et al. (2012) and Xue et al. (2015) proposed approaches to jointly model the probability of co-activation based on fMRI data while incorporating direct structural connections. They provide measures of functional co-activation deviating from standard measures of FC such as Pearson or partial correlations. Hinne et al. (2014) proposed a Bayesian approach which uses fMRI data to model the distribution of partial correlations for edges determined by the given SC information. Hence, functional connections only exist between anatomically connected regions which is an extremely restrictive assumption that does not reflect a realistic relationship between brain structure and function. Moreover, the above approaches use multi-subject data which requires registration of images to a shared template under the assumption that the volumes are similar and can be matched. Unfortunately, this assumption has limitations for human brain images considering the substantial variability in cortical anatomy and function (Zhu et al., 2012). This variability is especially pronounced during the developmental phases of childhood and adolescence as is the case with our motivating Philadelphia Neurodevelopmental Cohort (PNC) study. Ng et al. (2012) and Pineda-Pardo et al. (2014) proposed approaches for estimating sparse functional networks for individual subjects via an adaptive graphical lasso involving edge specific shrinkage parameters which are parametric functions of the SC strengths. Under these approaches, FC with less anatomical support are more heavily penalized, and vice-versa. However, the parametric form of the shrinkage parameters may not adequately capture the complex underlying structure-function relationships and does not account for heterogeneity in FC for a given SC strength resulting from non-anatomical sources of variation (see Figure 1 for a representative FC-SC relationship in the PNC study). Moreover, such a parametric relationship may lead to erroneous FC estimates when the SC information is mis-specified, rendering these approaches sensitive to the quality of anatomical information.

The above discussions highlight a serious need for integrative modeling approaches which adaptively estimate FC by incorporating structural knowledge in an appropriate manner. In designing such an approach, our primary goals for the method are that it will (a) provide accurate and reproducible brain network estimation, which is a topic of great importance in current literature (Varoquaux et al., 2010); and (b) specify a flexible structure-function relationship which is robust to mis-specification of SC information (arising from limitations in existing image acquisition technology) and can accommodate heterogeneity in FC for a given SC. We propose a novel hierarchical Bayesian Gaussian graphical modeling approach for estimating FC based on single subject fMRI data which incorporates given SC information in a manner that addresses the aforementioned objectives. The FC is computed via sparse precision matrices whose elements are estimated under Laplace type priors having edge specific shrinkage parameters that are random variables modeled using SC information and an independent baseline component. The prior encourages stronger FC given a large SC (and vice-versa), but also accounts for edge specific variations in FC unrelated to the brain anatomy, via the baseline component. The approach is thus flexible in accounting for anatomical knowledge with the FC being guided by, but not completely determined by, the SC information. Our method is motivated by the variable selection approach presented in Chang et al. (2016) which incorporates prior graph knowledge in a linear regression setting, but is distinct in addressing graphical model selection and precision matrix estimation, as well as the manner in which the prior knowledge is incorporated. Under certain choices of model parameters, the proposed approach reduces to an adaptive shrinkage approach specifying a parametric relationship between the shrinkage parameters and the anatomical information, similar to Ng et al. (2012) and Pineda-Pardo et al. (2014).

While Markov chain Monte Carlo (MCMC) can be used to implement the proposed approach, it is not scalable to large networks needed for whole brain connectome analysis as in our application. Moreover, MCMC samples cannot take exact zeroes under a Laplace prior, and an additional thresholding step is often needed for model selection. We propose an optimization algorithm to obtain the maximum a posteriori (MAP) estimate, which is computationally efficient, scales to a large dimensions, and does not require post-hoc thresholding. Under various simulation studies, we observe superior performance of the proposed method as compared to alternative approaches with or without SC information. The advantages of our approach become more evident as the degree of mis-specification of anatomical knowledge increases, and/or the number of nodes grows larger which is particularly relevant for whole brain connectome analysis.

Our efforts are motivated by data from the PNC study, a large-scale, NIMH funded initiative to understand the developmental trajectory of the brain from childhood to adolescence (Satterthwaite et al., 2014). The PNC data contains both DTI and resting state fMRI measurements from boys and girls ages 8-21 years, with suggestive but unclear SC-FC relationships (Figure 1). We fuse multimodal brain imaging data to examine gender differences in brain networks across three age brackets - pre-teens (ages 8-12), teens (ages 13-17), and young adults (ages 18-21) - and discover several gender based differences in FC within and across the age brackets. We also assess the reliability of computed network metrics across scanning sessions and find that the proposed approach yields strong reproducibility in the estimation of network metrics which is almost always higher than alternative approaches. While other studies have separately examined the reproducibility of functional and structural brain networks (see Welton et al. (2015) for a review), ours is one of the first to examine the reproducibility of anatomically informed functional networks to the best of our knowledge.

Section 2 describes the proposed methodology and the optimization routine for estimating networks, while sections 3 presents numerical studies and application of our method to PNC data. We conclude with a brief discussion in Section 4.

## 2 Materials and methods

### 2.1 Gaussian graphical model for brain networks

While early work on brain network estimation utilized Pearson correlation to measure undirected inter-regional dependencies, recent literature has focused on Gaussian Graphical Models (GGMs) which estimate a sparse inverse covariance matrix and measure functional connections via partial correlations (Smith et al., 2011). GGMs assume observations are normally distributed and that zeros in the inverse covariance matrix correspond to absent edges in the network, . A typical GGM specifies where represents the -dimensional pre-processed fMRI signal at time t across all ROIs, is the inverse covariance matrix (or the precision matrix) which belongs to the cone of symmetric positive definite matrices, denoted by . Under this framework, estimating is equivalent to estimating structural zeros in the positive definite precision matrix .

Due to tradeoffs between the cost and efficiency of information transfer, it is typically assumed that the brain seeks efficient organization favoring a sparse set of active connections at any point in time (Eavani et al., 2015). The GGM approach is well equipped to handle such sparse networks by imposing penalties that shrink sufficiently weak functional connections to zero, where the penalty under the graphical lasso (Friedman et al., 2008) is a popular choice (Ng et al., 2012; Monti et al., 2014; Pineda-Pardo et al., 2014). The graphical lasso can be thought of as an extension of the Lasso approach in regression settings and penalizes the full data likelihood to estimate the inverse covariance matrix as

(1) |

where is the sample covariance matrix, denotes the element-wise norm, and is the penalty parameter controlling overall network sparsity. If , one obtains the maximum likelihood estimate, while large values of shrinks an increasing number of off-diagonal elements to zero. The typical graphical lasso approach fits a series of graphs under various choices of the tuning parameter and chooses the optimal network as the one minimizing some goodness of fit criteria (Yuan and Lin, 2006).

### 2.2 Structurally informed Bayesian Gaussian graphical model

Bayesian GGM approaches have been successfully used for estimating brain networks (see Mumford and Ramsey (2014) for a review). One such approach is the Bayesian graphical lasso, which has similarities with the penalized graphical lasso approach in the sense that the maximum a posteriori (MAP) estimator is equivalent to the estimate under (1) (see Wang (2012)). The Bayesian graphical lasso, with a common shrinkage parameter, is defined as follows

(2) |

as in Wang et. al (2012), where the diagonal element is modeled under an exponential prior distribution , the off-diagonal element is modeled with double exponential or Laplace prior distribution , is the shrinkage parameter, and denotes the indicator function. In a fully Bayesian paradigm, is typically assigned a prior distribution, and is thus learned from the data, resulting in an adaptive shrinkage of the elements in .

In order to incorporate anatomical knowledge in functional connectivity estimation, we propose a hierarchical Bayesian structurally informed Gaussian graphical model (siGGM). It is based on the generic Bayesian GGM in (2), but involves edge specific shrinkage parameters which are modeled using anatomical knowledge. Throughout this article, we will denote the structural connectivity metric as for edge , where a larger value denotes a stronger anatomical connection and vice-versa. For example, in our data application, corresponds to the probability of SC obtained via probabilistic tractography. The proposed approach to estimating the brain functional network incorporating anatomical knowledge is defined as follows

(3) |

where denotes the prior distribution, denotes the collection of edge specific shrinkage parameters having a log-normal () type distribution, refers to the tuning parameter controlling the network’s overall sparsity and also corresponds to the scale parameter for the exponential prior on the diagonals, is a positive random variable which controls the average effect of SC on FC, denotes the random edge specific baseline component representative of non-anatomical sources of variations regulating FC, and is the intractable normalizing constant for the prior on the precision matrix depending on and . As in Wang et. al (2012), it can be shown that is finite and cancels out the normalizing constant in when computing the posterior distribution , resulting in a closed form expression which facilitates computation. We note that in the extreme case when , the model specifies a parametric relationship for fixed choices of and , which has similarities with Ng et al. (2012) and Pineda-Pardo et al. (2014).

The anatomically informed prior on the shrinkage parameters in (3) specifies a probabilistic relationship between the edge specific shrinkage parameters and the given SC knowledge via . In particular, increasing positive values of implies an increasing dependence on the given SC, potentially resulting in a functional connection even for small SC weights. Figure 2 illustrates that for large and increasing SC, the Laplace prior has heavier tails and less mass around zero, which i interpreted as increased probability of strong FC. In contrast, small values of do not result in a noticeable change in the prior distribution under varying SC strengths, implying a negligible relationship between SC and FC. Moreover, the shrinkage parameters are stochastically monotonically decreasing with respect to the SC strength, under the restriction . This implies that as the SC strength () for the edge is increased, the corresponding shrinkage parameter will take smaller values in probability, resulting in values of which are away from zero. This encourages the presence of FC at edge for large values of . Similarly small values of will encourage greater shrinkage for resulting in the absence of FC at edge .

Additionally, the baseline effect, , corresponds to variations in underlying neuronal activity that are independent of the brain anatomy. This formulation enables (a) more flexibility in the FC-SC relationship by allowing the possibility of strong FC when an anatomical connection is not obvious, and vice-versa; and (b) heterogeneity in FC across edges which possesses similar SC strength that is often encountered in practice. Overall, increasing(decreasing) absolute values of the baseline effect discourages(encourages) the presence of an edge in a manner that is independent of the anatomical information. The hyperparameters and are unknown and are learnt in a data-adaptive manner under the following priors

(4) |

where , and , are typically pre-specified. The choices of these hyperparameters are described in Appendix 5. We note that the scale parameter controls the overall network sparsity and is treated as a tuning parameter as in the penalized approach in (1), enabling the estimation of a series of networks with varying sparsity levels. The optimal network is chosen as the point estimate corresponding to the value of minimizing the Bayesian Information Criteria.

### 2.3 Model Estimation

Although the proposed model can be implemented using MCMC, it is not scalable to high dimensional settings involving a large number of nodes. Moreover, MCMC samples require a post hoc thresholding step to select important edges since estimates cannot take exact zeros under a Laplace prior. We bypass these limitations by computing a MAP estimate for the parameters of interest. Our iterative optimization approach employs an existing graphical lasso algorithm to sample the precision matrix given all other parameters, coupled with additional optimization steps to sample the shrinkage parameters and associated hyperparameters inherent in the Bayesian specification (3). In order to fit the proposed model, we estimate by maximizing the log-posterior distribution in (5), where and denotes the vector of edge-specific log-shrinkage parameters and baseline effects in (3), respectively. Note that the posterior distribution can be written as

We find the MAP solution for the model parameters by maximizing over the the posterior log-likelihood as , where

(5) | |||||

All parameters in the posterior distribution are updated iteratively until convergence. The precision matrix is updated given other parameters using the existing graphical lasso algorithm, whereas and are updated via closed form expressions and is updated via a Newton-Raphson step since a closed form solution does not exist. The iterative updates continue until for . At convergence, is the solution, where is the anatomically informed functional brain network based on single subject data. We note that one could alternatively treat and as tuning parameters and compute a range of networks over a grid of values, and then select the optimal network as one minimizing some goodness of fit criteria. However, this strategy did not result in adequate numerical performance, highlighting the advantages of specifying suitable priors on hyperparameters in order to estimate them in a data-adaptive manner. The computational steps for updating the model parameters are detailed in Appendix A.

## 3 Results and Discussion

### 3.1 Simulations

#### Simulation Setting

We conduct numerical studies to assess the performance of siGGM relative to SC naive and SC informed competitors. SC naive approaches are representative of methods that do not incorporate auxiliary information, and includes the graphical lasso or Glasso (Friedman et al., 2008), the partial correlation approach Space proposed by Peng et al. (2009), and the proposed approach in (3) without structural information obtained by setting , denoted siGGM(). SC informed approaches incorporate anatomical information into the estimation routine. We consider the Bayesian G-Wishart approach by Hinne et al. (2014) which treats FC as completely determined by SC and is denoted by G-Wishart and the adaptive graphical lasso approach by Pineda-Pardo et al. (2014) which specifies a parametric relationship between the shrinkage parameters and SC, and is denoted by aGlasso. All of the above approaches, except Space which estimates partial correlations, calculate sparse inverse precision matrices, where a zero off-diagonal entry implies the absence of an edge. The Glasso and Space approaches are implemented via the R packages glasso and space, respectively. We estimate the precision matrix under the G-Wishart approach using Matlab scripts available on the author’s website and incorporate adaptive weights in the glasso algorithm to implement aGlasso.

Data Generation: In order to assess the performance of our approach, we simulate data under three assumed network structures and consider various relationships between SC and FC. The network structures are (a) Erdos Renyi (ER) networks consisting of edges randomly generated with probability ; (b) small-world (SM) networks generated under the Watts-Strogatz model (Watts and Strogatz, 1998) in which most nodes may not be directly connected, but can reach other nodes via a small number of steps, and (c) scale-free (SF) networks generated using the preferential attachment model (Barabási and Albert, 1999), in which nodes are more likely to link to a highly connected node than to a node with few connections, resulting in a hub network. For each network, we consider varying the number of nodes corresponding to . The data was generated using a Gaussian graphical model with time points (), where is the precision matrix with zero off-diagonal elements corresponding to absent edges in the network . Once is determined based on the desired network structure, is constructed as follows: the non-zero off-diagonal elements corresponding to important edges are generated from a Uniform(-1,1) distribution and the diagonal elements were fixed to be one. In order to ensure positive definiteness, we subtracted the minimum of the eigenvalues from each diagonal element of the generated precision matrix.

Prior SC Knowledge: Conditional on , we construct several types of SC information, where the SC strengths were generated randomly between according to the following schemes. For scenario , we specify that 50% of those edges with strong FC also have strong SC(), while 25% of those edges with strong FC have moderate SC(), and the remaining 25% have weak SC(). For scenario , the proportion of edges that have strong FC, coupled with strong SC, moderate SC, and weak SC, are 30%, 35% and 35% respectively. For each of the scenarios and , we also consider two levels of mis-specification of the SC information, which are denoted as and respectively. For , we specify that 10% of those edges with zero FC have non-zero SC, while for we fix 20% of those edges with zero FC to have non-zero SC. All the other edges with zero or weak FC are assumed to have small SC, while remaining edges with moderate FC have non-zero SC strengths. We note that edges having zero FC but non-zero SC represent potential mis-specification of anatomical knowledge, based on the notion that strong SC typically underlies robust non-zero FC (Shen et al., 2015; Kemmer et al., 2017).

Comparison Metrics: To assess the performance under different approaches, we compute the area under the curve (AUC), which is a measure of the estimated sensitivity versus specificity over different network sparsity levels. Sensitivity is computed as , while specificity is defined as , where denote the number of edges that are true positives, true negatives, false positives and false negatives, respectively. To evaluate the point estimate of the network obtained under the Bayesian information criteria (BIC), we compute the Matthew’s correlation coefficient (MCC) which is a scalar measure combining sensitivity and specificity and is defined as MCC (Wang, 2012). We compute the relative norm error, , to assess accuracy in estimating the precision matrix encapsulating the FC strengths. Since brain networks are also often evaluated in terms of summary statistics reflecting network organization, we evaluate the accuracy in estimating the global efficiency which is a commonly used measure for global integration of brain connectivity.

#### Results

The results under siGGM and SC naive approaches under are presented in Table F.1, and Table F.2 displays results for SC informed approaches under various levels of mis-specification. Table F.1 illustrates that either the proposed siGGM approach, or the variant of the proposed approach with no prior knowledge (), have the lowest bias in estimating the global efficiency for all dimensions. Moreover, the proposed approach always has a higher MCC and AUC values, and a lower error norm compared to all other SC naive approaches. These results demonstrate the advantages of using structural knowledge to guide network estimation.

When the mis-specification levels are varied, Table F.2 illustrates that the proposed method has a consistently lower bias in estimating the global efficiency for both and mis-specification levels, compared to alternative SC informed approaches. Moreover, while the G-Wishart approach may have a higher MCC for when the mis-specification level is 10% (cases and ), the proposed method has a comparable or higher MCC for 20% mis-specification levels (cases and ). Moreover under , the MCC under the proposed approach is the highest for small-world and scale-free networks, and comparable to the G-Wishart method for the Erdos-Renyi network. We also note that while aGlasso often has the lowest MCC values, it may sometimes yield a higher AUC under small world and scale-free networks for . However, the proposed approach is shown to have the highest AUC for for all scenarios, highlighting the advantages of incorporating prior knowledge in a flexible manner in higher dimensions. Finally, siGGM consistently has the lowest error in estimating the precision matrix across all networks and dimensions. The above results illustrate a robust performance of the proposed method for and a superior performance for under the small world and scale-free networks, which closely resemble brain networks encountered in practical applications.

Although Table F.2 provides some idea about the relative performance under mis-specification, it is of interest to look at the effects of mis-specification in more detail. Hence, we examined the AUC and error values as the mis-specification level was gradually increased from 4% to 50%, under different networks for . The results, presented in Figure 4 illustrate that the proposed method consistently has higher AUC under the Erdos-Renyi network, and has a larger AUC for non-trivial mis-specification levels under the small world network. Further, the gains in AUC under the proposed approach seems to steadily increase with growing mis-specification levels under both these networks. Under the scale-free network, aGlasso has higher AUC under lower mis-specification levels when , but the proposed approach exhibits comparable or higher AUC as the mis-specification level is increased and it has a higher AUC for larger dimensions () as in Table F.2. For all networks, the proposed method is seen to have the lowest error across all mis-specification levels, while the error under the G-Wishart increases sharply as the mis-specification level is increased. These results provide clear evidence that while alternative SC informed approaches may occasionally perform better for lower mis-specification levels and moderate dimensions, the proposed siGGM method demonstrates a superior performance when the mis-specification level as well as the number of nodes is increased, which is of paramount importance in whole brain connectome studies involving high resolution atlases (Glasser et al., 2016; Power et al., 2011) and non-trivial mis-specification of the anatomical information.

Finally, we note that the siGGM can be implemented fairly quickly. On a 2.5Gz Intel Core i5 processor, the procedure estimates the optimal graph structure in approximately three seconds for p=40, twenty seconds for p=100, and approximately four minutes for p=200. While these computation times are slightly slower compared to generic graphical modeling approaches naive to anatomical knowledge, the overall computation is sufficiently quick and feasible for practical implementation in whole brain connectome analysis.

### 3.2 PNC Data Application

Existing literature has examined various neural substrates for age-related changes using structural and functional neuroimaging (Gur et al., 2012; Shaw et al., 2008; Raznahan et al., 2011). Moreover, gender differences have been extensively documented in behavioral measures (Halpern et al., 2007; Hines, 2010), structural neuroimaging (Lenroot et al., 2007), and functional imaging measures (Lenroot and Giedd, 2010). However, gender related differences in the developmental trajectory of the brain functional network from childhood to adolescence are still not understood well (Gur et al., 2012), and further, limited attempts have been made to investigate such differences by fusing functional and structural neuroimaging data. We use resting-state fMRI and DTI data from the Philadelphia Neurodevelopment Cohort (PNC) study to obtain preliminary answers to these questions. After estimating brain functional connectivity based on SC knowledge, we examine FC differences between boys and girls across different age groups.

We perform the analysis separately for each gender within the three age groups 8-12 (pre-teen), 13-17 (teen), 18-21 (young adult), where each age group contains approximately 9 to 12 individuals, and is constructed similarly to those in Ingalhalikar et al. (2014). All subjects are right-handed, physically, and mentally healthy, enabling a fair comparison between the groups. In addition to assessing gender based network differences, we also perform a secondary analysis to assess our method’s ability to reliably estimate functional networks. For this analysis, we split each subject’s resting-state fMRI time series into two equally sized scanning sessions (60 scans each) and calculate the intraclass correlation coefficient (refer to equation (5.3)) for seven network metrics which are widely used to summarize brain networks. These metrics include clustering coefficient, characteristic path length, local efficiency, global efficiency, modularity, hierarchy, and degree, and they were calculated with the Matlab toolboxes Brain Connectivity Toolbox (Rubinov and Sporns, 2010) and GRETNA (Wang et al., 2015).

#### Data preprocessing

Resting-state fMRI scans were acquired on a single-shot, interleaved multi-slice, gradient-echo, echo planar imaging (GE-EPI) sequence (Satterthwaite et al., 2014). Nominal voxel size is 3x3x3mm with full brain coverage achieved with the following parameters: TR/TE=3000/32 ms, flip=90, FOV=200 220 mm, matrix= 64 64, 46 slices, slice thickness/gap=3 mm/0 mm for a total of 6.2 minutes. Participants were instructed to remain awake, motionless, and fixated on a crosshair throughout the duration of the data acquisition. Several standard preprocessing steps were applied to the rs-fMRI data, including despiking, slice timing correction, motion correction, registration to MNI 2mm standard space, normalization to percent signal change, removal of linear trend, regressing out CSF, WM, and 6 movement parameters, bandpass filtering (0.009 to 0.08), and spatial smoothing with a 6mm FWHM Gaussian kernel. Subsequent voxel level data is aggregated into 90 regions of interest (ROI) based on the Automated Anatomical Labelling atlas (Tzourio-Mazoyer et al., 2002). For each ROI, the average time series of all constituent voxels represents the region’s temporal BOLD signal.

Diffusion weighted images permit us to localize and orient white matter fiber bundles via the diffusion of water in the brain. Images were acquired on a twice-refocused spin-echo (TRSE) single-shot EPI sequence for a total of 64 diffusion-weighted directions with b=1000 s/mm and 7 scans with b=0 s/mm (Satterthwaite et al., 2014). Acquisition parameters were TR/TE=8100/82ms, matrix=128128, FOV=240mm, slice thickness=2mm, GRAPPA factor=3. Due to gradient induced vibrations disturbing image quality, DWI images were acquired in two imaging runs to reduce the continuous duration in which subjects tolerate the scan. Standard pre-processing procedures, such as eddy current correction and bias-field correction are applied to the diffusion weighted data. Subsequently, we use the FSL functions bedpostx and probtracx2 to estimate the distribution of fiber tensors at each voxel and the count of white matter fibers tracts connecting all pairs of brain regions, respectively.

#### Results

In Figure 9 (A) and (B), we see that siGGM functional connectivity estimates are more highly correlated with non-zero anatomical connection strength across all subjects compared to Glasso and siGGM(), but have a lower correlation compared to the aGlasso approach specifying a parametric structure-function relationship. Additionally, we note in panels (C) and (D) of Figure 9 that siGGM leads to larger shrinkage and smaller variance for conditional dependencies between anatomically isolated brain regions compared to the generic graphical lasso without prior knowledge. This yields a smaller number of functional connections between anatomically disconnected ROIs. The above results indicate that our method adheres more closely to the brain’s anatomical connectivity without restricting the FC to be a deterministic function of structural pathways.

For network analysis, we classify each ROI into one of eight functional modules corresponding to resting state networks as defined in Smith et al. (2009). These functional modules include a medial visual network, an occipital pole and lateral visual network (“VIS”, 18 nodes), the default mode network (“DMN”, 8 nodes), a sensorimotor network (“SM”, 9 nodes), an auditory network (“AUD”, 10 nodes), an executive control network (“EC”, 19 nodes), right and left frontoparietal modules (“FPR” and “FPL”, 11 and 10 nodes, respectively) and an unknown module containing unassigned nodes (“UNK”, 5 nodes). Figure 12 shows that males and females have similar connectivity patterns with primarily positive connections within functional modules. Further comparisons of male and female brain networks within each age group reveals that consistent connections across age groups persist within module while inconsistent connections mainly exist between modules. After standardizing by the number of nodes in each module, the SM and AUD were found to be the two most highly connected functional modules in males and females across all age groups. Figure 15 (A) illustrates the similarity in network architecture for males and females with all metrics having non-significant differences across genders (except for local efficiency for teens and young adults), which implies shared patterns in brain organization across gender and age groups. Figure 15 (B) illustrates that males exhibit greater between module but smaller within-module connectivity differences in teens and young adults, which is supported by previous work and has been linked to variations in emotional identification and spatial cognitive tasks (Satterthwaite et al., 2016). However, we discover greater between- and within-module connections in pre-teen girls, which is an interesting finding that requires further examination.

As a second level of the analysis, we are also interested in the distribution of differentially weighted edges (DWE) between males and females, within each of the eight functional modules. DWE were identified as connections for which the FC strength was significantly different between genders under a permutation test. To evaluate if the number of DWEs within and between modules occur more often than allowed by chance, we define a goodness of fit measure (equation (5.2) in Appendix A) which represents the deviation between observed and expected numbers of DWE for each module block, standardized by the expected number. This measure captures whether a given module block has unusually high or low occurrence of DWE and enables us to identify modules with the most pronounced differences across gender. From the results presented in Table F.3, we discover statistically significant differences in the number of DWEs occurring in the executive control (EC) module in pre-teens and young adults, which is supported by previous results on gender related differences in the EC (Hyde, 1981; Mansouri et al., 2016). Table F.3 also suggests that gender based differences attenuate with development, with the largest number of DWEs in the pre-teen group and the smallest in the young adult group. We also find the DWE between the cingulum_ant_L in the EC and parietal_inf_L in the DMN exists in pre-teens, teens, and young adults, which suggests consistent gender based differences during the developmental phase. These regions are known to have brain volume differences between males and females which may point to subtle cognitive variations (Frederikse et al., 1999; Ruigrok et al., 2014).

A major challenge in resting state connectivity studies is to ensure reproducibility of the findings (Griffanti et al., 2016). We demonstrate that appropriately incorporating anatomical connectivity information leads to stable topological features of estimated networks across scanning sessions. Figure 16 displays the intraclass correlation coefficient (ICC) of seven network metrics under different approaches, where the details for computing the ICC are outlined in equation (5.3) in Appendix A. It is clear that the proposed siGGM produces estimates that have notably larger ICC measures for all the network metrics compared to all the other approaches considered. The reproducibility under the proposed approach is substantial for the clustering coefficient, global efficiency, and degree, and is moderate for all the other metrics. Moreover, it is reassuring to see that these three metrics with the highest ICC values under the proposed approach have been shown to be the most reproducible network metrics in independent studies (Niu et al., 2013; Telesford et al., 2010; Wang et al., 2011). In contrast, reproducibility is barely moderate under aGlasso for most metrics and weak under SC naive approaches. These findings highlight the benefits of incorporating anatomical information in a flexible manner. Although not presented, we note that the G-Wishart approach leads to an unrealistic ICC value of one in all cases, which is starkly different than the reliability values reported in previous studies (Welton et al., 2015). The perfect reliability is due to the fact that G-Wishart relies entirely on the SC information for specifying the functional network structure, resulting in the exact same network for both the sessions. Hence the reproducibility results under G-Wishart are not comparable.

## 4 Conclusion

In this paper, we introduce a novel Bayesian approach and an associated optimization algorithm for fusing structural and functional imaging data in order to estimate brain functional networks. We propose a flexible method for incorporating apriori known anatomical connectivity information in order to estimate the functional brain networks, bypassing the limitations of existing approaches by accommodating complex structure-function relationships while allowing unknown sources of variation independent of underlying anatomical structure. The proposed model is biologically more realistic compared to existing methods and often outperforms alternative approaches with or without anatomical knowledge as illustrated via extensive numerical studies. In particular, the advantages under the proposed method become more evident as the mis-specification levels for the anatomical knowledge and/or the number of nodes is increased, which has important practical implications. An analysis of the PNC data yields brain networks which have strongly reliable network metrics under the proposed approach, whereas the reproducibility under other approaches are moderate at best.

While we primarily focus on direct anatomical connections between brain regions, siGGM readily incorporates complex measures of anatomical connectivity corresponding to more indirect connections. In future work, we intend to explore various direct and indirect measures of structural connectivity so as to assess which anatomical measures yield the most meaningful and reproducible FC results.

## 5

### 5.1 Parameter Updates

We optimize (5) by iteratively updating model parameters as follows.

Update : Given the data, , and current estimates for all other model parameters, we solve

(5.1) |

for . This resembles the penalized likelihood framework of the traditional Gaussian graphical model. Define for and for . We can re-express (5.1) as

where we update using a quadratic approximation solver, QUIC, available in R.

Update : Given Y, , and , we update via the closed form equation

.

Update : Given Y and , we can update via closed form equation

where , , and

Update : Given Y, , , and , we can estimate for by solving

A closed form solution doesn’t exist, so we implement a Newton Raphson solver to find the optimal choice of . Re-expressing this problem, we have

where , denotes the upper diagonal elements of the structural connectivity matrix , is the element wise exponential for each component of , and 1 as a vector of 1’s of length . Since is symmetric and we do not shrink diagonal elements, we simplify our estimation of by only focusing upon the upper diagonal elements.

The Newton Raphson updating equation based on step size is where and , is a diagonal matrix with elements as the upper triangular elements of , and similarly for , and is an identity matrix. Since is a diagonal matrix, it is easily inverted and serves as an appropriate Hessian matrix. We search for the step size () using a back tracking line search for each update of as in Chang et. al (2017).

### 5.2 Hyper-parameter Choice and Initial Values

The proposed siGGM approach iteratively solves for the MAP estimator and works best when reasonable starting values are provided. We first find an initial estimate for the graph structure and the sparse inverse precision matrix (), using the graphical lasso. We initialize all edge specific penalty parameters as , which is the global tuning parameter corresponding to . We set corresponding to an uninformative prior which reflects our lack of knowledge regarding the baseline effects and choose as a default setting. We randomly generate the edge specific baseline effects from the prior distribution and use these as initial values. The initial value of is chosen by averaging , which is the average of all possible values under the relationship corresponding to . We choose .

Finally, we found that choosing and to attain E and Var allows incorporation of structural information in a flexible manner. However, larger first moments for the prior on may lead to increased false positives as our method places more weight on smaller structural connections, and similarly, smaller first moment may decrease the overall impact of structural information. For example, when and , we have E, which makes the siGGM indistinguishable from SC naive methods.

### 5.3 Measure for computing between module differences

We define the goodness of fit measure

(5.2) |

where are the indices corresponding to one of the functional modules, represents the observed number of DWEs in the block, represents the expected number of DWEs in the block when edges distribute randomly across the module blocks. measures the goodness of fit for each within-module block () or between-module block (). In equation (5.2), the expected value can be derived in a straightforward manner as for within module blocks () and for between-module blocks (), where represents the total number of nodes within the th module, and represents the proportion of DWE among all the edges across the network. Using 5000 permutations of group labels at each edge, the DWE are identified as those connections with significant FDR-adjusted p-values.

### 5.4 Calculation of ICC

The intraclass correlation coefficient is a widely used reliability metric for assessing test-retest reliability of brain network topology in neuroimaging applications. Using ICC(3,1), (two-way mixed single measures testing for consistency) we investigate the reliability of graph metrics across two scanning session (Guo et al., 2012; Telesford et al., 2010). The quantity is calculated as

(5.3) |

where is the number of scanning sessions per participant, is the between mean square and is the mean residual sum of squares. BMS captures the variability between subjects while EMS measures unexplained within-subject variation in functional connectivity across scanning sessions (see Shrout and Fleiss (1979)). This metric is commonly used to measure test-retest network stability in brain networks (Braun et al., 2012) with agreement scale (slight), (fair), (moderate), (strong), and (near perfect) as suggested by Telesford et al. (2010).

## 6

Pre-Teen | ||||||||
---|---|---|---|---|---|---|---|---|

Unknown | Visual | DMN | SM | Aud | EC | FP left | FP right | |

Unknown | 0 | |||||||

Visual | 8 | 24 | ||||||

DMN | 7 | 15 | 0 | |||||

SM | 3 | 14 | 3 | 2 | ||||

Aud | 2 | 17 | 13 | 5 | 4 | |||

EC | 6 | 32 | 9 | 20 | 21 | 42 | ||

FP left | 2 | 14 | 8 | 4 | 10 | 13 | 0 | |

FP right | 5 | 18 | 9 | 11 | 11 | 16 | 3 | 6 |

Teen | ||||||||

Unknown | Visual | DMN | SM | Aud | EC | FP left | FP right | |

Unknown | 2 | |||||||

Visual | 6 | 22 | ||||||

DMN | 3 | 8 | 2 | |||||

SM | 1 | 11 | 9 | 4 | ||||

Aud | 3 | 7 | 7 | 13 | 4 | |||

EC | 8 | 22 | 13 | 12 | 12 | 18 | ||

FP left | 3 | 16 | 7 | 6 | 9 | 13 | 0 | |

FP right | 2 | 12 | 9 | 5 | 15 | 18 | 8 | 2 |

Young Adult | ||||||||

Unknown | Visual | DMN | SM | Aud | EC | FP left | FP right | |

Unknown | 2 | |||||||

Visual | 1 | 10 | ||||||

DMN | 1 | 14 | 2 | |||||

SM | 1 | 10 | 9 | 4 | ||||

Aud | 1 | 15 | 5 | 5 | 8 | |||

EC | 5 | 16 | 10 | 12 | 12 | 26 | ||

FP left | 1 | 7 | 5 | 7 | 11 | 9 | 6 | |

FP right | 1 | 8 | 7 | 9 | 7 | 17 | 8 | 0 |

M.A.(2014), ‘Structurally-informed bayesian functional connectivity analysis’, NeuroImage 86, pp. 294-305.

### Footnotes

- Research reported in this publication was supported by the National Institute Of Mental Health of the National Institutes of Health under Award Number ROI MH105561 and R01MH079448. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

### References

- Barabási, A.-L. and Albert, R. (1999), ‘Emergence of scaling in random networks’, science 286(5439), 509–512.
- Bressler, S. L. and Tognoli, E. (2006), ‘Operational principles of neurocognitive networks’, International Journal of Psychophysiology 60(2), 139–148.
- Chang, C., Kundu, S. and Long, Q. (2016), ‘Scalable bayesian variable selection for structured high-dimensional data’, arXiv preprint arXiv:1604.07264 .
- Damoiseaux, J. S. and Greicius, M. D. (2009), ‘Greater than the sum of its parts: a review of studies combining structural connectivity and resting-state functional connectivity’, Brain Structure and Function 213(6), 525–533.
- Eavani, H., Satterthwaite, T. D., Filipovych, R., Gur, R. E., Gur, R. C. and Davatzikos, C. (2015), ‘Identifying sparse connectivity patterns in the brain using resting-state fmri’, Neuroimage 105, 286–299.
- Frederikse, M. E., Lu, A., Aylward, E., Barta, P. and Pearlson, G. (1999), ‘Sex differences in the inferior parietal lobule’, Cerebral Cortex 9(8), 896–901.
- Friedman, J., Hastie, T. and Tibshirani, R. (2008), ‘Sparse inverse covariance estimation with the graphical lasso’, Biostatistics 9(3), 432–441.
- Glasser, M. F., Coalson, T. S., Robinson, E. C., Hacker, C. D., Harwell, J., Yacoub, E., Ugurbil, K., Andersson, J., Beckmann, C. F., Jenkinson, M. et al. (2016), ‘A multi-modal parcellation of human cerebral cortex’, Nature 536(7615), 171–178.
- Griffanti, L., Rolinski, M., Szewczyk-Krolikowski, K., Menke, R. A., Filippini, N., Zamboni, G., Jenkinson, M., Hu, M. T. and Mackay, C. E. (2016), ‘Challenges in the reproducibility of clinical studies with resting state fmri: An example in early parkinson’s disease’, Neuroimage 124, 704–713.
- Guo, C. C., Kurth, F., Zhou, J., Mayer, E. A., Eickhoff, S. B., Kramer, J. H. and Seeley, W. W. (2012), ‘One-year test–retest reliability of intrinsic connectivity network fmri in older adults’, Neuroimage 61(4), 1471–1483.
- Gur, R. C., Richard, J., Calkins, M. E., Chiavacci, R., Hansen, J. A., Bilker, W. B., Loughead, J., Connolly, J. J., Qiu, H., Mentch, F. D. et al. (2012), ‘Age group and sex differences in performance on a computerized neurocognitive battery in children age 8- 21.’, Neuropsychology 26(2), 251.
- Halpern, D. F., Benbow, C. P., Geary, D. C., Gur, R. C., Hyde, J. S. and Gernsbacher, M. A. (2007), ‘The science of sex differences in science and mathematics’, Psychological science in the public interest 8(1), 1–51.
- Hermundstad, A. M., Bassett, D. S., Brown, K. S., Aminoff, E. M., Clewett, D., Freeman, S., Frithsen, A., Johnson, A., Tipper, C. M., Miller, M. B. et al. (2013), ‘Structural foundations of resting-state and task-based functional connectivity in the human brain’, Proceedings of the National Academy of Sciences 110(15), 6169–6174.
- Hines, M. (2010), ‘Sex-related variation in human behavior and the brain’, Trends in cognitive sciences 14(10), 448–456.
- Hinne, M., Ambrogioni, L., Janssen, R. J., Heskes, T. and van Gerven, M. A. (2014), ‘Structurally-informed bayesian functional connectivity analysis’, NeuroImage 86, 294–305.
- Honey, C. J., Thivierge, J.-P. and Sporns, O. (2010), ‘Can structure predict function in the human brain?’, Neuroimage 52(3), 766–776.
- Hyde, J. S. (1981), ‘How large are cognitive gender differences? a meta-analysis using! w and d..’, American Psychologist 36(8), 892.
- Ingalhalikar, M., Smith, A., Parker, D., Satterthwaite, T. D., Elliott, M. A., Ruparel, K., Hakonarson, H., Gur, R. E., Gur, R. C. and Verma, R. (2014), ‘Sex differences in the structural connectome of the human brain’, Proceedings of the National Academy of Sciences 111(2), 823–828.
- Kemmer, P. B., Bowman, F. D., Mayberg, H. and Guo, Y. (2017), ‘Quantifying the strength of structural connectivity underlying functional brain networks’, arXiv preprint arXiv:1703.04056 .
- Lenroot, R. K. and Giedd, J. N. (2010), ‘Sex differences in the adolescent brain’, Brain and cognition 72(1), 46–55.
- Lenroot, R. K., Gogtay, N., Greenstein, D. K., Wells, E. M., Wallace, G. L., Clasen, L. S., Blumenthal, J. D., Lerch, J., Zijdenbos, A. P., Evans, A. C. et al. (2007), ‘Sexual dimorphism of brain developmental trajectories during childhood and adolescence’, Neuroimage 36(4), 1065–1073.
- Mansouri, F. A., Fehring, D. J., Gaillard, A., Jaberzadeh, S. and Parkington, H. (2016), ‘Sex dependency of inhibitory control functions’, Biology of sex differences 7(1), 11.
- Messé, A., Rudrauf, D., Benali, H. and Marrelec, G. (2014), ‘Relating structure and function in the human brain: relative contributions of anatomy, stationary dynamics, and non-stationarities’, PLoS computational biology 10(3), e1003530.
- Monti, R. P., Hellyer, P., Sharp, D., Leech, R., Anagnostopoulos, C. and Montana, G. (2014), ‘Estimating time-varying brain connectivity networks from functional mri time series’, NeuroImage 103, 427–443.
- Mumford, J. A. and Ramsey, J. D. (2014), ‘Bayesian networks for fmri: a primer’, Neuroimage 86, 573–582.
- Ng, B., Varoquaux, G., Poline, J.-B. and Thirion, B. (2012), ‘A novel sparse graphical approach for multimodal brain connectivity inference’, Medical Image Computing and Computer-Assisted Intervention–MICCAI 2012 pp. 707–714.
- Niu, H., Li, Z., Liao, X., Wang, J., Zhao, T., Shu, N., Zhao, X. and He, Y. (2013), ‘Test-retest reliability of graph metrics in functional brain networks: a resting-state fnirs study’, PLoS One 8(9), e72425.
- Peng, J., Wang, P., Zhou, N. and Zhu, J. (2009), ‘Partial correlation estimation by joint sparse regression models’, Journal of the American Statistical Association 104(486), 735–746.
- Pineda-Pardo, J. A., Bruña, R., Woolrich, M., Marcos, A., Nobre, A. C., Maestú, F. and Vidaurre, D. (2014), ‘Guiding functional connectivity estimation by structural connectivity in meg: an application to discrimination of conditions of mild cognitive impairment’, Neuroimage 101, 765–777.
- Power, J. D., Cohen, A. L., Nelson, S. M., Wig, G. S., Barnes, K. A., Church, J. A., Vogel, A. C., Laumann, T. O., Miezin, F. M., Schlaggar, B. L. et al. (2011), ‘Functional network organization of the human brain’, Neuron 72(4), 665–678.
- Raznahan, A., Lerch, J. P., Lee, N., Greenstein, D., Wallace, G. L., Stockman, M., Clasen, L., Shaw, P. W. and Giedd, J. N. (2011), ‘Patterns of coordinated anatomical change in human cortical development: a longitudinal neuroimaging study of maturational coupling’, Neuron 72(5), 873–884.
- Rubinov, M. and Sporns, O. (2010), ‘Complex network measures of brain connectivity: uses and interpretations’, Neuroimage 52(3), 1059–1069.
- Ruigrok, A. N., Salimi-Khorshidi, G., Lai, M.-C., Baron-Cohen, S., Lombardo, M. V., Tait, R. J. and Suckling, J. (2014), ‘A meta-analysis of sex differences in human brain structure’, Neuroscience & Biobehavioral Reviews 39, 34–50.
- Satterthwaite, T. D., Connolly, J. J., Ruparel, K., Calkins, M. E., Jackson, C., Elliott, M. A., Roalf, D. R., Hopson, R., Prabhakaran, K., Behr, M. et al. (2016), ‘The philadelphia neurodevelopmental cohort: a publicly available resource for the study of normal and abnormal brain development in youth’, Neuroimage 124, 1115–1119.
- Satterthwaite, T. D., Elliott, M. A., Ruparel, K., Loughead, J., Prabhakaran, K., Calkins, M. E., Hopson, R., Jackson, C., Keefe, J., Riley, M. et al. (2014), ‘Neuroimaging of the philadelphia neurodevelopmental cohort’, Neuroimage 86, 544–553.
- Shaw, P., Kabani, N. J., Lerch, J. P., Eckstrand, K., Lenroot, R., Gogtay, N., Greenstein, D., Clasen, L., Evans, A., Rapoport, J. L. et al. (2008), ‘Neurodevelopmental trajectories of the human cerebral cortex’, Journal of Neuroscience 28(14), 3586–3594.
- Shen, K., Mišić, B., Cipollini, B. N., Bezgin, G., Buschkuehl, M., Hutchison, R. M., Jaeggi, S. M., Kross, E., Peltier, S. J., Everling, S. et al. (2015), ‘Stable long-range interhemispheric coordination is supported by direct anatomical projections’, Proceedings of the National Academy of Sciences 112(20), 6473–6478.
- Shrout, P. E. and Fleiss, J. L. (1979), ‘Intraclass correlations: uses in assessing rater reliability.’, Psychological bulletin 86(2), 420.
- Smith, S. M., Fox, P. T., Miller, K. L., Glahn, D. C., Fox, P. M., Mackay, C. E., Filippini, N., Watkins, K. E., Toro, R., Laird, A. R. et al. (2009), ‘Correspondence of the brain’s functional architecture during activation and rest’, Proceedings of the National Academy of Sciences 106(31), 13040–13045.
- Smith, S. M., Miller, K. L., Salimi-Khorshidi, G., Webster, M., Beckmann, C. F., Nichols, T. E., Ramsey, J. D. and Woolrich, M. W. (2011), ‘Network modelling methods for fmri’, Neuroimage 54(2), 875–891.
- Sporns, O. (2013), ‘Structure and function of complex brain networks’, Dialogues in clinical neuroscience 15(3), 247.
- Telesford, Q. K., Morgan, A. R., Hayasaka, S., Simpson, S. L., Barret, W., Kraft, R. A., Mozolic, J. L. and Laurienti, P. J. (2010), ‘Reproducibility of graph metrics in fmri networks’, Frontiers in neuroinformatics 4.
- Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B. and Joliot, M. (2002), ‘Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain’, Neuroimage 15(1), 273–289.
- Varoquaux, G., Gramfort, A., Poline, J.-B. and Thirion, B. (2010), Brain covariance selection: better individual functional connectivity models using population prior, in ‘Advances in neural information processing systems’, pp. 2334–2342.
- Venkataraman, A., Rathi, Y., Kubicki, M., Westin, C.-F. and Golland, P. (2012), ‘Joint modeling of anatomical and functional connectivity for population studies’, IEEE transactions on medical imaging 31(2), 164–182.
- Wang, H. (2012), ‘Bayesian graphical lasso models and efficient posterior computation’, Bayesian Analysis 7(4), 867–886.
- Wang, J.-H., Zuo, X.-N., Gohel, S., Milham, M. P., Biswal, B. B. and He, Y. (2011), ‘Graph theoretical analysis of functional brain networks: test-retest evaluation on short-and long-term resting-state functional mri data’, PloS one 6(7), e21976.
- Wang, J., Wang, X., Xia, M., Liao, X., Evans, A. and He, Y. (2015), ‘Gretna: a graph theoretical network analysis toolbox for imaging connectomics’, Frontiers in human neuroscience 9.
- Watts, D. J. and Strogatz, S. H. (1998), ‘Collective dynamics of "small-world" networks’, nature 393(6684), 440–442.
- Welton, T., Kent, D. A., Auer, D. P. and Dineen, R. A. (2015), ‘Reproducibility of graph-theoretic brain network metrics: a systematic review’, Brain connectivity 5(4), 193–202.
- Xue, W., Bowman, F. D., Pileggi, A. V. and Mayer, A. R. (2015), ‘A multimodal approach for determining brain networks by jointly modeling functional and structural connectivity’, Frontiers in computational neuroscience 9.
- Yuan, M. and Lin, Y. (2006), ‘Model selection and estimation in regression with grouped variables’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
- Zhu, D., Li, K., Guo, L., Jiang, X., Zhang, T., Zhang, D., Chen, H., Deng, F., Faraco, C., Jin, C. et al. (2012), ‘Dicccol: dense individualized and common connectivity-based cortical landmarks’, Cerebral cortex 23(4), 786–800.