Small Worldness in Dense and Weighted Connectomes

Colon-Perez, Luis M.; Couret, Michelle; Triplett, William; Price, Catherine C.; Mareci, Thomas H.

doi:10.3389/fphy.2016.00014

Original Research ARTICLE

Front. Phys., 10 May 2016 | https://doi.org/10.3389/fphy.2016.00014

Small Worldness in Dense and Weighted Connectomes

Luis M. Colon-Perez¹^*,

Michelle Couret²,

William Triplett³,

Catherine C. Price⁴ and

Thomas H. Mareci³

¹Department of Psychiatry, University of Florida, Gainesville, FL, USA
²Department of Medicine, Columbia University in the City of New York, New York, NY, USA
³Department Biochemistry and Molecular Biology, University of Florida, Gainesville, FL, USA
⁴Department Clinical and Health Psychology, University of Florida, Gainesville, FL, USA

The human brain is a heterogeneous network of connected functional regions; however, most brain network studies assume that all brain connections can be described in a framework of binary connections. The brain is a complex structure of white matter tracts connected by a wide range of tract sizes, which suggests a broad range of connection strengths. Therefore, the assumption that the connections are binary yields an incomplete picture of the brain. Various thresholding methods have been used to remove spurious connections and reduce the graph density in binary networks. But these thresholds are arbitrary and make problematic the comparison of networks created at different thresholds. The heterogeneity of connection strengths can be represented in graph theory by applying weights to the network edges. Using our recently introduced edge weight parameter, we estimated the topological brain network organization using a complimentary weighted connectivity framework to the traditional framework of a binary network. To examine the reproducibility of brain networks in a controlled condition, we studied the topological network organization of a single healthy individual by acquiring 10 repeated diffusion-weighted magnetic resonance image datasets, over a 1-month period on the same scanner, and analyzing these networks with deterministic tractography. We applied a threshold to both the binary and weighted networks and determined that the extra degree of freedom that comes with the framework of weighting network connectivity provides a robust result as any threshold level. The proposed weighted connectivity framework provides a stable result and is able to demonstrate the small world property of brain networks in situations where the binary framework is inadequate and unable to demonstrate this network property.

Introduction

The brain is a heterogeneous system comprised of a broad range of white matter (WM) connection strengths [1, 2]. Computational models of neuronal networks account for this heterogeneity by weighting the network edges. Chavez et al. found that weighted edges enhanced synchronization between network units [3]. To account for this enhanced synchronization and the broad range of large and small white matter tracts in the brain, the analysis of brain networks must include a connectivity weighting parameter to reflect the real underlying organizational structure of the brain [4, 5]. In addition to producing more realistic network models of the brain, weighting the edges adds an extra degree of freedom in calculated network parameters.

The topological organization of structural and functional brain networks has been extensively studied in vivo with magnetic resonance imaging (MRI) [6–10] and graph theory provides an appropriate framework in which to elucidate the topological organization of brain networks [4]. This approach has revealed several basic network characteristics of the brain, such as high clustering, short path lengths, modularity, and small world organization. In order to generate brain structural networks, white matter structural connections (network edges) are estimated from diffusion-weighting MRI (dMRI) using tractography between gray matter regions (network nodes). The resulting networks suggest only a small world organization for brain networks [11–14]. But these studies generally employ binary connectivity representations, which assume every connection between nodes is equivalent [11, 15], so the connection either exists or does not exist, because the strength of connection is not used to differentiate between connections. Weighted networks can differentiate the strength of connections, however weighted connectome studies have focused on the difference of specific metrics between healthy controls and pathological subjects. Weighted connectome studies have shown increased changes in node strength, efficiency, clustering between patients with various neurological diseases compared to healthy controls [16–20]. These studies illustrate that global changes may take place in the organization of network connections in neurological disorders. This suggests that the small-worldness topological index may also show these organizational changes in weighted connectomes.

Complex network models have been employed to describe various real world networks [21, 22] and weighted network parameters have been introduced, but a comprehensive weighted connectivity framework has not been reported which allows the estimation of brain network topological features, such as small worldness [23]. The usual binary framework approach to estimating the topology of brain networks includes the calculation of network metrics; degree distribution, path length and clustering coefficients [23]. In separate published works, all of these metrics have been generalized to their weighted counterparts. In the next section of this paper, these generalized metrics are reviewed then used to estimate the topological features (small worldness) of brain networks. This weighted network approach leads to a more realistic model of the brain network and a more robust characterization of the brain network topology.

Network Metrics

The adjacency matrix, A, is the most basic mathematical representation of network connectivity [24]. For simple and undirected graphs, the adjacency matrix is square and symmetric (i.e., N × N, where N is the number of nodes in the network):

\begin{array}{l} A = (\begin{matrix} a_{11} & \dots & a_{1 n} \\ ⋮ & ⋱ & ⋮ \\ a_{n 1} & \dots & a_{n n} \end{matrix}), & (1) \end{array}

Unweighted networks, known as binary networks, are described in an adjacency matrix by either the presence or absence of an edge connecting any two nodes where these binary edges have the value,

\begin{matrix} a_{i j} = {\begin{array}{l} 1, i f a n e d g e c o n n e c t s n o d e s n_{i} a n d n_{j} \\ 0, \end{array} . & (2) \end{matrix}

In a weighted networks, the adjacency matrix is modified by adding a weighting parameter to represent the connection strength of each edge,

\begin{array}{l} a_{i j}^{'} = a_{i j} \cdot w (e_{i j}) . & (3) \end{array}

Here the edge weight, w(e_ij), characterizes the connectivity strength between each pair of nodes. The weighted adjacency matrix allows the calculation of network metrics, analogous to those derived from the binary adjacency matrix [25], and is the foundation of a framework to study the topology of weighted brain networks. Since translational diffusion is antipodal symmetric, only undirected graphs will be considered for dMRI-derived connections; hence, the elements of the binary and the weighted adjacency matrix will be symmetric (i.e., a_ij = a_ji). A simple overall measure of network connectivity, the graph density (ρ), is defined as the ratio of the number of edges in the network graph to the maximum number of possible edges. This parameter is used to quantify gross alterations in networks resulting from changes, such as a selected threshold.

Node Connectivity

The node degree, k_i, for the i^th node, n_i, represents the number of nodes connected to node, n_i, and is calculated from the binary adjacency matrix as follows,

\begin{array}{l} k_{i} = \sum_{j = 1}^{N} a_{i j} . & (4) \end{array}

A degree distribution, expressed as the probability that nodes have a degree, k, in the network of N nodes can be used to provide a simple representation of connectivity in the network. In weighted networks, a parameter analogous to degree (Equation 4) is the node connection strength, s_i [26],

\begin{array}{l} s_{i} = \sum_{j = 1}^{N} a_{i j}^{'} = \sum_{j = 1}^{N} a_{i j} \cdot w (e_{i j}) . & (5) \end{array}

Like the degree distribution, the node connection strength distribution provides a simple representation of connectivity in a weighted network.

Mean Geodesic Path Length

Another important network metric is the geodesic path length. Geodesic paths are useful because they represent efficient pathways for the transmission of information within the network. In binary networks, a geodesic path length, d_ij, is the smallest number (dimensionless) of edges required to connect node i to node j [25]. The mean geodesic path length for a specific node in the brain has been associated with the efficiency of the overall network structure [27]. For an undirected binary network, the mean geodesic path length, ${\hat{d}}_{i}$ , for node, n_i, is defined as,

\begin{array}{l} {\hat{d}}_{i} = \frac{\sum_{j = 1, i \neq j}^{N} d_{i j}}{N - 1} . & (6) \end{array}

Starting from the definition of the binary geodesic path length, a geodesic path length for weighted networks, $d_{i j}^{w}$ , can be defined using a Dijkstra's algorithm [28]: first define the paths with the same and smallest number of edges between nodes n_i and n_j, then choose the path with highest sum of edge weights. An illustration of the algorithm is shown in the graph of Figure 1. Using the weighted adjacency matrix, the algorithm starts by determining the low-cost path between nodes connected by one edge. For the path between nodes n₁ and n₂, the original form of Dijkstra's algorithm would assign the lowest-cost (highest connection strength) to the path along the edges from n₁ to n₃, and n₃to n₂ because the total connection strength is stronger than the direct edge between node n₁ and n₂. For brain networks, Dijkstra's algorithm is modified to give highest priority to the connections with the fewest number of edges then sort those to find the lowest cost. Therefore, the low-cost path between nodes n₁ and n₂ is the direct edge between node n₁ and n₂. As the modified Dijkstra's algorithm continues, the low-cost path is determined for nodes connected by two or more edges. For example, the low-cost path from node n₁ to n₆ is found by starting at node n₁, which has three neighboring nodes, and two of these nodes have a two-edge path, which connect to n₆. The algorithm compares $d_{16}^{w}$ for the path w(e₁₂) + w(e₂₆) and w(e₁₃) + w(e₃₆). With the edge weight represented by the thickness of the edge, a visual inspection of each these paths in Figure 1 shows that w(e₁₂) + w(e₂₆) < w(e₁₃) + w(e₃₆); therefore in this example, the low-cost path is $d_{16}^{w} = w$ (e₁₃) + w(e₃₆). In general, the process of determining the low-cost paths starting at node n_i begins with directly connected nodes and is repeated for path lengths with increasing number of steps until all of the (N − 1) low-cost paths are determined. This results in the following form of the mean (low-cost) geodesic weighted-path length for node i;

\begin{array}{l} {\hat{d}}_{i}^{w} = \frac{\sum_{(j = 1, i \neq j)}^{N} d_{i j}^{w}}{N - 1} . & (7) \end{array}

FIGURE 1

Figure 1. Eleven-node weighted network. Edge thickness represents the relative edge weight strength, e.g., e₁₃ > e₁₄.

In the context of brain networks, a step along the mean geodesic weighted-path is associated with the strongest connection (i.e., highest sum of edge weights) between nodes.

Clustering Coefficient

The motifs of connections within the network provide information about network structure and how the network functions. The clustering coefficient, c_i, is a measure of the motif of triangular connections for the n_i node, which quantifies the local connectivity around this node by measuring the connectivity between the neighboring nodes directly connected to node i. The clustering coefficient in binary form is given by [23, 24],

\begin{array}{l} c_{i} = \frac{2 E_{j m}}{k_{i} (k_{i} - 1)} = \frac{2}{k_{i} (k_{i} - 1)} [\frac{1}{2} \sum_{j, m = 1}^{N} a_{i j} a_{j m} a_{m i}] \\ = \frac{1}{k_{i} (k_{i} - 1)} \sum_{j, m = 1}^{N} a_{i j} a_{j m} a_{m i}, & (8) \end{array}

where E_jm (defined in square brackets) is the number of edges connecting the neighbors of node i (i.e., the number of connected closed triangles, where the three nodes involving node i are fully connected). The normalization k_i(k_i − 1)/2 represents the number of possible connections among the neighbors of node i. In the brain, the clustering coefficient has been associated with specialized processing (e.g., sensory input analysis, such as visual and auditory) [27] where nearby nodes work together to achieve complex tasks. Therefore, high clustering of the neighboring nodes ultimately allows efficient communication and complex task processing.

Weighted network generalizations of the binary clustering coefficient in Equation (8) have been reviewed by Saramaki et al. [29], who conclude that none of the proposed generalizations provide an all-purpose weighted clustering coefficient, and they suggest the characterization of a network using a clustering coefficient should be made from two perspectives: binary and weighted. Onnela et al. [30] introduced a weighted clustering coefficient, c_{i, O} in Equation (9) below, that replaces the sum of triangular binary adjacency coefficients with the scaled edge weights (scaled by the maximum edge weight in the network, as shown in Equation 11 below).

\begin{matrix} c_{i, O} = \frac{1}{k_{i} (k_{i} - 1)} \sum_{j, m = 1}^{N} [\hat{w} (e_{i j}) \hat{w} (e_{j m}) \hat{w} (e_{m i})]^{1 / 3}, & (9) \end{matrix}

Also Zhang and Horvarth [31], introduced a weighted clustering coefficient, c_{i, Z} in Equation (10) below, using the sum over the product of each triangle of normalized edge weights in the numerator, but included a normalization factor in the dominator using scaled edge weights.

\begin{array}{l} c_{i, Z} = \frac{\sum_{j, m = 1}^{N} ŵ (e_{i j}) ŵ (e_{j m}) ŵ (e_{m i})}{{(\sum_{j = 1}^{N} ŵ (e_{i j}))}^{2} - \sum_{j = 1}^{N} ŵ {(e_{i j})}^{2}}, & (10) \end{array}

where

\begin{array}{l} ŵ (e_{i j}) = \frac{w (e_{i j})}{max (w (e))} . & (11) \end{array}

As ŵ(e_ij) → 1 (or w(e_ij) → 1) in all triangles, both of these weighted clustering coefficients approach c_i, thus ensuring that the results in weighted networks will yield the expected binary results when the weighted edges are converted to binary edges. In this study, these weighted clustering coefficients will be compared to the binary clustering coefficient, as suggested by Saramaki et al. [29].

Small Worldness

Many real-world networks have been shown to have high overall clustering and short average path lengths, which reflect regional specialization and information transfer efficiency [32], and are classified as small world networks. Small world networks are more clustered than random networks, defined with the Erdös' and Rényi's model of random graphs [33], yet display similar geodesic path lengths. Humphries and Gurney [34] introduced a small-worldness metric which uses the network-average clustering coefficient and path length values, relative to these metrics for random networks, to provide an overview of connectivity in the entire network. The small-worldness metric is obtained by taking the ratio of two ratios: (1) the ratio of the average clustering coefficient (Equations 8–10) of the brain network, c_g, to clustering coefficient, c_N, of a random networks,

\begin{array}{l} γ = \frac{c_{g}}{c_{N}}, & (12) \end{array}

and (Equation 2) the ratio of average geodesic path length [average of all nodes for the path of Equations (6 and 7)] of brain networks, l_g, to the average geodesic path length of a random networks, l_N,

\begin{array}{l} λ = \frac{l_{g}}{l_{N}} . & (13) \end{array}

Then the small-worldness metric, sw, is given by the ratio of Equations (12 and 13),

\begin{array}{l} s w = \frac{γ}{λ} . & (14) \end{array}

In a binary framework, small world networks have γ > 1 and λ ~ 1 and present a high degree of connections among the neighbors of any node, compared to a random network, but small average path lengths are preserved. In weighted networks, any node may present a strong level of connectivity among its neighbors, while preserving small path lengths, so that the small-worldness metric can be used to characterize weighted networks.

Materials and Methods

MRI Acquisition

The University of Florida Institutional Review Board approved this human study. One healthy subject was scanned 10 times over the course of 1 month, which provided 10 controlled brain data set from which to determine the network properties across different MRI acquisitions [35, 36]. In all 10 sessions the subject was awake and scans were acquired roughly at the same time of the day (~7:00 a.m.). The subject was scanned in a 3 T Siemens Verio system in the Shands Hospital of the University of Florida [37]. High angular resolution diffusion imaging (HARDI) data was obtained with a spin echo preparation and an echo planar imaging [38] readout using the following parameters: TR/TE = 17300/81 ms, two scans without diffusion weighting, six diffusion gradient directions with b-value of 100 s/mm², and 64 diffusion gradient directions with b-value of 1000 s/mm². The diffusion gradients were distributed following a scheme of electrostatic repulsion [39]. The diffusion-weighted images covered the entire brain with an isotropic resolution of 2.0 mm, field of view (FOV) of 256 × 256 mm and 73 slices. To increase the spatial resolution, data were interpolated with cubic convolution [40] using the CONGRID function in the Interactive Data Language (IDL; Exelis Visual Information Systems, McLean, VA) to a voxel dimension of 1 mm³. In addition, a high-resolution T1 structural scan of the entire brain was acquired with TR/TE = 2500/3.77 ms, resolution 1 mm isotropic, FOV of 256 × 256 mm and 176 slices.

Post Processing

The diffusion weighted scans were corrected for eddy current distortion using FSL's eddy_correct algorithm [41]. FA and AD maps were created with in-house software using IDL. The displacement probability function (DPF) was calculated for each voxel using a mixture of Wishart distribution (MOW) [42] implemented with in-house, C-based software. The MOW distributions were obtained using a distribution of positive definite matrices (rank-2 tensors) that allow the estimation of multiple fiber orientations in each MR voxel. Then the fiber orientations in each voxel are identified as the DPF maxima with probability of 0.5 or larger.

Using the MOW fiber orientations, deterministic tractography was performed to generate streamline maps throughout the entire brain. Tractography was performed following a modified version of the fiber assignment by continuous tracking (FACT) algorithm [43], in which the direction with the least angular deviation along the incoming fiber path was selected to be continued at each iteration of the tractography process. Fiber maps were calculated using an in-house, C-based software. Seed points were placed uniformly throughout the MR voxels to avoid the ambiguity of random placement of seed points. From each seed point, one streamline is launched bi-directionally for each estimated displacement probability maxima contained in the voxel. Each streamline is propagated by stepping half of the MR voxel width in the direction that is most in line with the streamline's current direction of propagation. Tractography was performed with 125 seeds per voxel, a fiber step size of 0.5 of the voxel size, using the following criteria: no step-to-step track deviations greater than 50°, and no tracking into voxels with low anisotropy values (FA < 0.05). If the estimated track does not meet these criteria, the streamline is stopped.

Network Construction

Brain networks were generated from 68 anatomical nodes segmented [44, 45] using the structural T1-weighted image with an automatic segmentation algorithm in FreeSurfer (Laboratory for Computational Neuroimaging, A. A. Martinos Center for Biomedical Imaging, Charlestown, MA). The structural image and diffusion weighted images was spatially registered using FSL's FLIRT [46] algorithm, using an affine transformation. Employing FLIRT's transformation matrix output, the FreeSurfer derived nodes were then registered to DWI using FSL's ApplyXFM.

The network edges are defined by streamlines connecting any two nodes and the edge weight is defined as [37],

\begin{array}{l} w (e_{i j}) = (\frac{V_{v o x e l}}{P_{v o x e l}}) (\frac{2}{A_{i} + A_{j}}) \sum_{p = 1}^{P_{v o x e l}} \sum_{m = 1}^{M} \frac{χ_{R} (f_{p, m})}{l (f_{p . m})}, & (15) \end{array}

where

\begin{matrix} χ_{R} (f_{p, m}) = {\begin{matrix} 1, f_{p, m} \in R \\ 0, f_{p, m} \notin R \end{matrix}, & (16) \end{matrix}

V_voxel is the MR voxel volume, P_voxel is the number of seed points per voxel,A is the surface area of each node, M is the number of voxels making up the edge, f_{p, m} is the streamline originating from seed point p in voxel m, l(f_{p, m}) is length of f_{p, m} streamline making up the edge, and χ_R(f_{p, m}), is the characteristic function for the set R of streamlines connecting nodes n_iand n_j. The characteristic function removes (i.e., filters) streamlines that originate within the nodes and voxels that do not represent the WM path connecting the nodes [37]. This form of the edge weight calculation eliminates the effect of length of the streamlines, seeding paradigm, and tractography-specific experimental factors from the calculation [37]. The inverse sum along with the characteristic function (Equation 16) yields the number of streamline fibers connecting the nodes, and the normalization by the average surface area of the nodes reduces the bias introduced by node size. This form of the edge weight allows the quantification of white matter connection strength between any two nodes in a dimensionless and scale invariant manner.

In brain networks created from tractography, thresholds are usually employed to reduce the number of false edges and the effects of voxel volume averaging in network measures. In binary networks, thresholding to produce a sparse networks (i.e., the number of edges is much less than 30% of the total number of possible edges) is used to estimate better the topological properties of binary networks, since dense networks are highly connected and may lead to an almost complete graph (every node connected to all other nodes in the graph). But thresholding binary networks introduce two problems in networks calculations. Firstly the threshold selection is arbitrary and secondly networks with different thresholds are not independent samples [47]. Arbitrarily chosen thresholds do not ensure that small and real edges are included or that artifactual edges are completely eliminated, since there is no clear criterion to achieve the elimination of artifactual streamlines. Also, the statistical results obtained for all binary networks will depend on the selected thresholds, since binary networks created at different thresholds are pseudo replicas of group-level networks [47]. This means that the network results obtained from the same data, expressed as a sparse network, are also obtained in a dense network, but not vice-versa. This indicates that there is direct dependence of the calculated network properties on the applied threshold, which leads to different results at different thresholds. In this paper, a threshold parameter is selected that can be applied to both binary and weighted networks. Then network metrics (discussed in Section Network Metrics) are calculated at different threshold levels in order to determine whether weighted networks allow better estimation of an appropriate threshold or if it is possible to eliminate thresholds from the analysis.

In this work we will use the streamline fiber count is a threshold parameter. For the weighted network discussed here, a high seed density (P_voxel = 125) is used to maximize the accuracy of the calculated edge weights [37]. But weighted-network thresholds should be set low enough so that “weak” edges are not removed just to create sparse networks. To illustrate how the streamline-fiber-count threshold approach can be used, assume that the “smallest” edge that will contribute to the brain network has a length of 2 mm and streamline fibers that are entirely contained within two voxels. Using parameters from Equation (15), a maximum of 250 streamlines (M × P = 2 × 125 = 250) should be connecting the nodes. Therefore, a threshold of 25 streamlines would retain all edges that have at least ~10% (i.e., 25/250) of the streamlines that can potentially make this “smallest” edge. Therefore, using the number of streamlines connecting any two nodes as a threshold parameter, and not the edge weight value, will allow this threshold to be applied to edges of both the binary and weighted brain network. Considering the high volume of seeds-per-voxel used in this study, four threshold values will be used to examine the effects of thresholding on calculated binary and weighted network parameters. First the network metrics will be calculated without applying any threshold, i.e., any two nodes connected by at least one streamline will have an edge, then edges that contain at least 25, 50, and 125 or more streamlines will be considered valid edges.

All network metrics described in the previous section are calculated with in-house code written in R (http://cran.us.r-project.org/) and the network package in R (http://cran.r-project.org/web/packages/network/index.html) was used to create and modify networks, and create relational data within the R interface. Three dimensional brain networks were displayed with BrainNet [48].

Null Hypothesis Graphs

Network metrics are influenced by the decisions made when the network is generated (e.g., identification of nodes and edges). Therefore, to test the significance of differences in the metrics (i.e., clustering coefficients and path lengths) between binary and weighted brain networks, these metrics are compared to the same metrics obtained for a null hypothesis network [49]. A null-hypothesis binary-network is formed by assigning edges at random using two constraints: (1) the number of nodes and (2) the degree distribution remains identical to the original brain network. Since both the brain network and null hypothesis network are constrained by the number of nodes and the degree distribution, differences between them are not due to local differences in connectivity (i.e., node degree discrepancies) [7].

To create weighted null hypothesis networks, a third constraint is used to account for the additional degree of freedom resulting from weighting the edges. This third constraint requires that the edge weight distributions in the null hypothesis and original brain network are the same [14]. The process of creating the null hypothesis network starts by storing lists of the edge weight values and the node degree distribution (binary) for each brain network [49]. Then starting with the same number of nodes as the brain network, the edges are randomly connected in node pairs, while preserving the node degree distribution. As each edge is assigned, an edge weight is randomly selected from the edge weight list to create the null hypothesis network. Ultimately, this creates a random network with the same node degree distribution (i.e., binary null hypothesis network) and edge weight distribution (i.e., weighted null hypothesis network) as the calculated brain network.

Results

Both binary and weighted brain networks are shown in Figure 2 at the threshold of 125 or more fibers. The mean binary graph density at different thresholds for all 10 networks, derived from the 10 dMRI measurements of a single subject, are shown in Figure 3. As the threshold increases, the graph density is reduced from ~51% (0 threshold) to ~30% (125 threshold). The level of variation (error bars in Figure 2) is consistent across all thresholds with coefficients of variation from 3.5 to 3.8%. From this point on, the order of all threshold comparisons will be from 0, 25, 50, to 125.

FIGURE 2

Figure 2. Brain networks. Representative (A) binary and (B) weighted network. The nodes in the binary network are scaled by the node degree and in the weighted network by the connection node strength.

FIGURE 3

Figure 3. Graph density vs. threshold. Density is obtained by taking the ratio of the number of edges in a graph to the total possible number of edges in the graph. The threshold is the number of steamline fibers below which an edge is removed from the network.