Minimum spanning tree analysis of the human connectome.

One of the challenges of brain network analysis is to directly compare network organization between subjects, irrespective of the number or strength of connections. In this study, we used minimum spanning tree (MST; a unique, acyclic subnetwork with a fixed number of connections) analysis to characterize the human brain network to create an empirical reference network. Such a reference network could be used as a null model of connections that form the backbone structure of the human brain. We analyzed the MST in three diffusion-weighted imaging datasets of healthy adults. The MST of the group mean connectivity matrix was used as the empirical null-model. The MST of individual subjects matched this reference MST for a mean 58%-88% of connections, depending on the analysis pipeline. Hub nodes in the MST matched with previously reported locations of hub regions, including the so-called rich club nodes (a subset of high-degree, highly interconnected nodes). Although most brain network studies have focused primarily on cortical connections, cortical-subcortical connections were consistently present in the MST across subjects. Brain network efficiency was higher when these connections were included in the analysis, suggesting that these tracts may be utilized as the major neural communication routes. Finally, we confirmed that MST characteristics index the effects of brain aging. We conclude that the MST provides an elegant and straightforward approach to analyze structural brain networks, and to test network topological features of individual subjects in comparison to empirical null models.

INTRODUCTION

The human brain is a complex network and its structural organization is closely related to its function. Network analysis has increased our understanding of the brain, including the genetic basis of connection formation, the effects of maturation and aging on neural network organization, and the correlates of high level cognitive processing (Bullmore and Sporns, 2009; Bullmore and Sporns, 2012; Fornito, Zalesky, & Bullmore, 2016). In network studies, the brain is analyzed as a set of nodes (gray matter regions) and their connections. The connections can be structural, such as white matter tracts, functional, for example, representing synchronized activity between populations of neurons, or theoretical, based on correlations in characteristics such as cortical thickness or gene expression. One of the key characteristics of human brain networks is the presence of so‐called hub nodes, which are central brain regions with a relative high number of connections (Gong et al., 2009; van den Heuvel and Sporns, 2011). A general backbone of cortico‐cortical highways has been described in the healthy human brain using graph theoretical analysis of diffusion weighted imaging tractography, which includes connections between a subset of hub nodes, the so‐called rich club (Hagmann et al., 2008; van den Heuvel, Kahn, Goni, & Sporns, 2012). Mounting evidence suggests that structural brain networks, particularly their backbone, are altered in neurological and psychiatric disorders, and heterogeneity in these alterations is evident between diseases (Bullmore and Sporns, 2009; Crossley et al., 2014; Stam, 2014).

Although our knowledge of brain networks is rapidly growing, several methodological issues remain unresolved. One of these issues is the lack of a reference network, which describes the topology of a “standard” human brain; such a reference network is lacking even if scanning and preprocessing conditions are fixed (Fornito, Zalesky, & Breakspear, 2013; Stam et al., 2014). A reference could be useful to compare results of different studies and of different populations, provided that the preprocessing pipeline, node and edge definitions are similar. A reference network could also serve as a null‐model to characterize network alterations in individual subjects, and to develop theoretical models of altered connectivity patterns during maturation or due to disease (Fornito et al., 2013). Although an “average” brain network has been described, the use of this network as a reference to study differences in specific conditions has been hindered by another major issue in brain network studies, namely the so‐called “thresholding problem.” The thresholding problem refers to the large effect of trivial characteristics of the network on conventional graph measures, such as the number and strength of the connections (van Wijk, Stam, & Daffertshofer, 2010). The problem can be illustrated as follows: let's consider two networks A and B, where the number of connections in network A is twice the number of connections in network B. We now wish to compare the efficiency of both networks, using the shortest path length between the two most remote nodes in the network. In this comparison, network A is likely to be more efficient than network B in an absolute sense, simply because there are more connections in network A. However, if we were to normalize or correct for the greater number of connections in network A, the efficiency of both networks might be comparable. Therefore, an absolute difference in efficiency does not necessarily imply that the topological organization of the connections is also more efficient in one of the networks (Stam et al., 2014; van Wijk et al., 2010). Proportional thresholding is frequently applied to weighted connectivity matrices to overcome this problem and yield networks that are matched in terms of total number of connections; network characteristics are then based on this thresholded network containing a subset of strong connections. However, this approach affects the graph properties in an unpredictable manner, because the impact of setting a threshold varies, depending on the underlying network topology (Stam et al., 2014; van den Heuvel et al., 2017; van Wijk et al., 2010). Handling of the thresholding problem is likely to be one of the reasons that empirical studies have reported contradictory findings, for example of both increased and decreased characteristic shortest path length in patients with Alzheimer's Disease (Fornito et al., 2013; Tijms et al., 2013). Decreased shortest path lengths in the presence of an overall reduction in connectivity strength in a patient group may be due to the inclusion of spurious long‐distance connections in patients due to the use of a more lenient proportional threshold.

We propose to study the backbone structure of the brain network by mapping a mathematically defined subset of connections known as the minimum spanning tree (MST; Figure 1; Kruskal, 1956; Wang, Hernandez, & Van Mieghem, 2008). The MST is a subnetwork containing the strongest connections from the set of all available weighted connections. It connects all the nodes in the network without forming cycles or loops. The MST is unaffected by the thresholding problem provided that the ranking of edge weights remains unaltered, while its characteristics can be interpreted along the lines of conventional graph measures (Stam et al., 2014; Tewarie, van Dellen, Hillebrand, & Stam, 2015). Importantly, the MST is a unique subset of connections provided that all connection weights are unique, which makes it possible to create an empirical reference network of the connections in a healthy human brain. Recent studies have shown that MST analysis can be used to capture network alterations due to aging and disease in both functional and structural imaging data (Boersma et al., 2010; Otte et al., 2015; van Dellen et al., 2015, 2014).

Figure 1

The concept of the minimum spanning tree. Three minimum spanning tree network types. (a) shows a path tree, where every node except the two end nodes or leafs (red) is connected to its two neighbors (low leaf number), but it takes a lot of steps to reach the other end of the network (high diameter). (c) shows a star tree, which consists of a central node (high betweenness centrality, in green) that is connected to all other nodes (high degree), which are all leaf nodes. This network is highly efficient (low diameter), but may result in an overload of information flow through the central hub node. (b) represents a hierarchical tree, which is a possible intermediate between the two extremes [Color figure can be viewed at http://wileyonlinelibrary.com]

The topology of a network depends not only on the definition of connections, but also on the definition of nodes, in this case the brain regions that are considered to be part of the connectome (Zalesky et al., 2010). It remains unclear how brain network analysis may be influenced by the use of atlases to define cortical regions, or the exclusion of subcortical brain regions.

In the present work, we studied whether MST analysis could be used to define a reference or standard backbone of the human connectome. Structural brain networks were reconstructed in three datasets of healthy adults using diffusion weighted imaging both including and excluding subcortical regions. We defined a reference MST based on the average weighted connectivity matrix of each population of healthy subjects, and we tested the variance in overlap with this reference network for individual subjects. We studied the effect of using different scanners at 1.5 T and 3 T, the effect of different scanning sites, and processing pipelines. The topology of a network depends not only on the definition of connections, but also on the definition of nodes (Zalesky et al., 2010). We therefore studied the effects of using different brain atlases for node definitions, and the exclusion of subcortical brain regions. To investigate biological relevance in healthy adults, we aimed to replicate previous findings on the relation between MST topology and age and gender. Finally, we tested whether the MST captured key aspects of structural network organization of the human brain, including the consistent presence of core connections across subjects, and the location of hub nodes and rich club nodes.

METHODS

Dataset 1 (Netherlands dataset)

Subjects

Data of 46 healthy subjects were analyzed. The population consisted of 19 males and 27 females (41% and 59%, respectively), with a mean age of 39 years (Standard Deviation (SD) 14 years, range 20–65). Thirty‐three (72%) subjects were right‐handed and 13 (28%) were non‐right‐handed. Subjects were recruited as control subjects for a study in the department of psychiatry at the University Medical Center Utrecht via a website (http://www.verkenuwgeest.nl.) Inclusion criteria were: (a) no current or past psychiatric disorders; (b) no chronic somatic disorder; (c) absence of alcohol or drug abuse for at least 3 months. See de Weijer et al. (2013) for a detailed description of the inclusion procedure.

Image acquisition

MRI scans were acquired on a 3 T Philips Achieva with an eight‐channel SENSE head‐coil as previously reported in (de Weijer et al., 2013). Two sets of Diffusion Weighted Imaging scans and a T1‐weighted scan for anatomical reference were collected from each participant. Two sets of transverse DWI scans were obtained to increase the signal‐to‐noise ratio, using the following parameters: a single shot EPI‐DTI scan consisting of 30 diffusion weighted scans (b = 1000s/mm2) with noncolinear gradient directions and an average of five diffusion unweighted scans (b = 0 s mm2), TR/TE = 7,035/68 ms, FOV = 240 mm, matrix 128 × 128, 75 slices with 2 mm thickness, no gap, SENSE factor = 3, EPI factor = 35. The second set was identical to the first but acquired with reversed k‐space readout (anterior direction) to allow correction of geometric EPI distortions in the image‐processing step. The diffusion‐weighted scans were used for reconstruction of the fiber tracts. The parameters for the T1‐weighted scan were: TR/TE = 9.87/4.6 ms, flip angle = 8°, FOV 224 × 160 × 168, matrix = 256 × 256, slice thickness 1 mm (no gap).

Image processing

The T1‐weighted image was used for anatomical reference and for network node definition. Brain regions (network nodes) were selected automatically using the FreeSurfer segmentation pipeline (V5.1; http://surfer.nmr.mgh.harvard.edu) (Fischl et al., 2004). The brain was divided into 82 distinct anatomical regions of the Desikan‐Killiany Atlas using an automated segmentation procedure, consisting of 14 subcortical structures, and a parcellation of the neocortex into 68 regions. An individual mask was created containing all 82 regions for each subject. This mask was registered to the diffusion tensor weighted data set using a six‐parameter rigid body transformation with nearest‐neighbor interpolation.

Preprocessing of the diffusion‐weighted scans was performed with the diffusion toolbox of Andersson and Skare (2002) and Andersson, Skare, and Ashburner (2003) using an in‐house developed software pipeline (Mandl et al., 2010). A distortion map was calculated based on five averaged b = 0 images and applied to the two sets of 30 direction‐weighted images to correct for susceptibility artifacts. This resulted in a corrected DWI set consisting of a single b = 0 image and 30 corrected weighted images, thereby avoiding the need for nonlinear registration approaches to the T1‐weighted images (Andersson et al., 2003). The DWI set was corrected for eddy‐current distortions and small head movements by realigning all scans to the diffusion‐unweighted image (Andersson and Skare, 2002). The diffusion tensors where obtained using M‐estimators to limit the influence of possible outliers (Chang, Jones, & Pierpaoli, 2005).

Tracts were reconstructed in native space. All possible tracts in the brain were reconstructed individually using the diffusion weighted images with an in‐house implementation of the fiber assignment by continuous tracking (FACT) algorithm (Mandl et al., 2010) with the following parameter settings: 8 seed points per voxel, stopping criteria: minimum voxel fractional anisotropy = 0.1 or tract angle between streamline steps > 45°, and maximum average angle with neighboring voxels = 45°. After tractography was completed, all fibers shorter than 10 mm were discarded, as they were deemed spurious. Finally, all remaining fibers were linearly extended by 5 mm in the orientation prior to termination to maximize the probability of penetration into the grey matter. The presence of a white matter connection between two grey matter regions was determined by labeling each streamline with the grey matter areas it touched based on the anatomical segmentation mask. The number of streamlines (NOS) was summed between all possible node pairs and represented in an 82 × 82 structural connectivity matrix M nos.

Network reconstruction

Further analyses were performed using Matlab (The Mathworks,Inc., Natick, MA, United States), and the brain connectivity toolbox (https://sites.google.com/site/bctnet). Networks were reconstructed for cortical and subcortical ROIs (N = 82), and for cortical ROIs only (N = 68). Several methods for connection estimation based tractography data have been proposed (Jones and Leemans, 2011). We used the number of streamlines to calculate the structural connectivity matrix M1. However, network reconstruction based on the number of streamlines for segmented regions may be affected by the size of a region of interest (ROI), as regions with a bigger surface may have a higher probability of being adjacent to a streamline (Fornito et al., 2013; van den Heuvel and Sporns, 2011). We therefore also computed a structural connectivity matrix, where the number of streamlines was weighted by regional volume (M volume_weighted_nos), and one where the FA values were used as edge weights M FA (van den Heuvel and Sporns, 2011).

To quantify how MST characteristics may be influenced by connection definitions, we calculated correlations between nodal MST properties and ROI volume and surface, and estimated the MST edge length based on the Euclidean distance between the centroid voxel of each ROI.

Dataset 2 (Australia dataset)

Participant data were obtained from the Australian Schizophrenia Research Bank (ASRB), which is an Australian resource and storage facility of medical research data collected across 5 Australian states and territories. Exclusion criteria included any organic brain disorder, history of brain trauma followed by a long period of amnesia (>24 h), mental retardation (full‐scale IQ < 70), movement disorders, current drug or alcohol dependence, a personal or family history of psychosis or bipolar I disorder, and electroconvulsive therapy in the past 6 months.

Participants were English speaking, predominately of European ancestry and aged between 18 and 65 years (Loughland et al., 2010). For this study, participants with both a diffusion‐weighted imaging (DWI) and structural magnetic resonance imaging scan (sMRI) scan that was processed using Freesurfer were selected. This comprised 197 healthy controls (99 male). All participants provided informed consent for the analysis of their stored data. Study procedures were approved by the Melbourne Health Human Research Ethics Committee. Data use was granted through a variation of ASRB Access Agreement P49 (Zalesky).

T1‐weighted magnetic resonance imaging scans (sMRI) and diffusion‐weighted imaging (DWI) scans were acquired in each participant with a Siemens Avanto 1.5‐T system (Siemens, Erlangen, Germany) across five different sites in Australia. The same acquisition sequences for T1‐weighted and DWI scans was used across all sites. Full details related to sequence parameters have been described elsewhere (Klauser et al., 2017).

The imaging data were preprocessed with the FMRIB Software Library (FSL; Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012). In brief, preprocessing involved eddy correction for motion and gradient‐induced currents using affine registration to the first volume. The resulting transformation was inverted and applied to the automated anatomical atlas (AAL; Tzourio‐Mazoyer et al., 2002).

Automated whole‐brain fiber tracking was performed in native diffusion space for each individual to generate streamlines tracing out the trajectories of all major cortico‐cortical and cortico‐subcortical fiber bundles. For each subject, five thousand streamlines were seeded from randomly positioned coordinates comprising each brain region of the AAL‐90 atlas. Each streamline was propagated using the second‐order integration over fiber orientation distributions (iFOD2) probabilistic tractography algorithm (Tournier, Calamante, & Connelly, 2010). Streamlines were guided by fiber orientations inferred using spherical deconvolution with a maximum harmonic order (l max) of 4 (Tournier, Calamante, & Connelly, 2007). Propagation was terminated if either a minimum angle threshold of 45° was violated or the streamline propagated beyond the white matter mask. The strength of a connection (connectivity) between a pair of brain regions was defined by the NOS. These steps were performed using the MRtrix package (Tournier, Calamante, & Connelly, 2012) and Matlab.

Validation connectivity matrix

We compared the overlap of cortical connections in the MST of the Netherlands dataset to a reference MST based on a publicly available connectivity matrix from the International Consortion for Brain Imaging (ICBM) dataset (http://umcd.humanconnectomeproject.org). All DTI images in this sample were downloaded from the ICBM database (Mazziotta et al., 2001).

Data consisted of 138 healthy adults (72 female, 66 male), age range 18–90 years old, (mean 43.7, SD 15.2), obtained using a 1.5 T Siemens system, 30 diffusion directions, b = 1000 s/mm2, voxel size = 2.5 mm3 DTI images were preprocessed for eddy correction and tensor estimation using the FSL Diffusion Toolbox (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT). The probability distribution of fiber directions in each voxel was estimated using BEDPOSTX (Behrens, Berg, Jbabdi, Rushworth, & Woolrich, 2007) allowing two crossing fibers within each voxel. The dyads for the first and second vectors of diffusion direction within each voxel were used for as the input for deterministic tractography, using the fiber assignment by continuous tracking (FACT) algorithm in Diffusion Toolkit (http://trackvis.org/dtk). Whole brain tractography was carried out, propagating fibers from each voxel with a maximum turn angle of 50°. Fibers were smoothed using a spline filter and all fibers < 5 mm were excluded. MP‐RAGE scans had a voxel size = 1 mm3. Five repeat scans per subject were input into the Freesurfer pipeline where they were parcellated using the recon‐all program into 68 cortical regions of the Desikan‐Killiany Atlas (Desikan et al., 2006) and their corresponding subcortical white matter counterparts. The MP‐RAGE scan was registered to the eddy‐corrected DTI B0 image with FSL FLIRT (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT) using an affine transformation with 7 degrees of freedom. The resultant transformation matrix was used to register the 68 gray and white cortical parcellations to DTI space. We then counted number of fibers that originated in one given region and terminated in another, repeating this process between all pairs of regions. These raw fiber counts were then scaled by the mean volume of each pair of regions in order to obtain the fiber density count, the value used to populate the subject's connectivity matrix. Finally, these 68 × 68 matrices were averaged across subjects to obtain the group‐level structural connectivity matrix. The MST was calculated from this group average connectivity matrix and will be referred to as MSTICBM.

Minimum spanning tree analysis

The connectivity matrices were analyzed as weighted graphs using NOS as edge weights. Each region of the atlas was considered a node and all connections between any pair of ROIs were considered as edges. We used MST analysis to characterize network topology from the connectivity matrices using Kruskal's algorithm (Kruskal, 1956). This procedure starts with ranking all connection weights from lowest weight to highest weight. Since we are interested in the strongest connections, we ranked all connections from highest to lowest weight (formally, our procedure therefore reconstructs the maximum spanning tree). We start by disconnecting all nodes, and add the connection with the highest weight. Next, the connection with the second highest weight is added and this procedure is repeated until all nodes are connected. If adding a new connection results in a cycle or loop, this connection is discarded, and the next connection ranked by weight is selected.

The resulting MST is a subgraph of the underlying network, which connects all nodes such that the connections with minimum cost are included, but without forming loops (Wang et al., 2008). The MST is unique, provided that the connection weights are also unique. We applied this procedure to reconstruct the tree containing the connections with the maximum number of fibers, resulting in an MST for each subject containing a fixed number of 82 nodes and 81 links. In this way, it is guaranteed that no differences between connection density or strength are present between subjects. The full analysis pipeline is visualized in Figure 2.

Figure 2

Analysis pipeline. Data analysis pipeline. The brain was parcellated in cortical and subcortical regions of interest (ROIs). Connections were estimated with DTI tractography and connection strength was based on the number of streamlines or fractional anisotropy. The minimum spanning tree of each subject was reconstructed from this structural connectivity matrix. In addition, a mean connectivity matrix of all subjects was calculated, and the minimum spanning tree of this connectivity matrix was used as a reference network. Networks of individual subjects were then compared to this reference matrix [Color figure can be viewed at http://wileyonlinelibrary.com]

The weighted connectivity matrix was transformed into an unweighted MST (i.e., binary graph containing edge weights of 0 and 1 only; all MST edges were set to 1), which was used to further characterize network topology. A reference MST (MSTref) was calculated based on the average connectivity matrix of all subjects. In this average matrix, connections were given the streamline count averaged over all subjects. The similarity between the MST of each individual subject (MSTsubject) and MSTref was quantified by calculating the fraction of edges that was present in both MSTs (MSToverlap; ranging between 0 (no matching edges) and 1 (exact match) (Lee, Kim, & Jung, 2006). The use of an average matrix for construction of the MSTref may potentially be biased by outliers when connections in the MST are absent in only a few subjects (de Reus and van den Heuvel, 2013). We therefore also calculated the occurrence matrix of all possible connections in the MST of individual subjects (sum MSTsubject). For a pair of regions, the occurrence matrix indicates the number of subjects with an MST that has a connection between those regions. The MST of this occurrence matrix was computed, the MSToccurrence, and the overlap between MSToccurrence and MSTref was used to estimate the effect of outliers. Several characteristics of the MST of each subject and of MSTref were calculated to quantify network topology: diameter, leaf fraction, tree hierarchy and kappa, and maximum nodal degree and betweenness centrality (Boersma et al., 2012; Tewarie et al., 2015). Nodal characteristics betweenness centrality and degree were analyzed to describe the role of specific ROIs in the network. A detailed description and definitions of these measures are given in Table 1.

Table 1

Concepts and terminology

Characteristic	Definition	Interpretation	Formula
Degree	Number of links for a given node	Measure of regional importance. Nodes with a high degree may be considered “hubs,” i.e., crucial regions on the functional brain network	k_i=∑j∈Na_ij N is the number of nodes
Betweenness centrality (BC)	BC of a node u is defined as the number of shortest paths between any two nodes i and j in the network that are passing u, divided by the total number of shortest paths. BC ranges between 0 (leaf node) and 1 (central node in a star‐like network). The BC of the tree was characterized by the BCmax, i.e., the BC of the node with the highest BC in the tree.	Nodes with a high BC are considered “hub nodes” not based on their number of connections, but on their importance for global communication in the network. Maximum BC describes the importance of the most central node, which is a measure of central network organization.	BC_i=1(n−1)(n−2)∑h,j∈Nh≠j, h≠i, ≠iρhj(i)ρ_hj ρhj is the number of shortest paths between h and j, and ρhj(i) is the number of shortest paths between h and j that pass through i.
Diameter	Characterizes the largest distance between any two nodes, normalized for the total number of connections: D = d/M.	Measure of the efficiency of global network organization. In a network with a low diameter, information is efficiently processed between remote brain regions.	D=d/M M is the total number of links or maximum possible leaf number
Leaffraction	Measure based on the leaf number, which is defined as the number of nodes that have only one connection. It ranges between 2 (a line‐topology; such a tree is called a path) and a maximum value M = n − 1 (with n the number of nodes) (a star‐like topology). Leaf fraction is the leaf number divided by the maximum possible leaf number: L f = L/M.	Measure of global network topology that describes to what extent the network has a central organization. When the leaf fraction is high, communication is largely dependent on hub nodes.	L_f=L/M
Degree divergence (κ)	Measure of the broadness of the degree distribution.	Related to resilience against attacks of complex networks. Higher values of kappa reflect a broader degree distribution, and higher vulnerability for targeted attacks.	κ=〈k^2〉〈k〉
Tree hierarchy	Characterizes a hypothesized optimal topology of efficient organization while preventing information overload of central nodes	For a line‐like topology T h = 0, for a star‐like topology T h = 0.5, and for trees with a configuration between these 2 extreme situations, T h can have values of T h → 1.	T_h=12mBC_max
Overlap	The fraction of links that two MSTs (MSTx and MSTy) have in common. This value can range between 0 and 1.		Overlap=MST_x∩MST_yM

MST measures and their definitions (Boersma et al., 2012; Tewarie et al., 2014b).

ρ is the number of shortest paths between h and j, and ρ(i) is the number of shortest paths between h and j that pass through i.

The overlap between hub nodes and rich club nodes was analyzed. Hub nodes were defined as nodes with an MST degree or MST betweenness centrality of at least 1 standard deviation above the mean for the MST, or nodes with a degree of M + 1 SD in the full, unweighted connectivity matrix (corresponding to a degree k > 13 in our dataset) (van den Heuvel and Sporns 2011). We established the presence of rich club architecture in our dataset for individual subjects with a degree ranging from k = 14 to k = 18. Rich club analyses were performed on the unweighted network; edges were binarized and were considered present for edge weights > 0. The rich club coefficient ϕ of the full, unweighted matrix of each individual subject was calculated and normalized by comparison to ϕrandom, which was defined as the average rich club coefficient of 100 randomized networks with preserved degree distribution (van den Heuvel and Sporns 2011). We found a mean ϕnormalized = 1.05 for k = 14; ϕnormalized = 1.13 for k = 18; and ϕnormalized = 1.05–1.13 for k 15–17; rich club architecture is considered to be present for ϕ > 1. We followed van den Heuvel and Sporns in our definition of rich club nodes: a group mean unweighted connectivity matrix was calculated and edges were included when present in >75% of the subjects. Rich club nodes were nodes with a degree of mean + 1 SD (van den Heuvel & Sporns, 2011).

Statistical analysis

Statistical analyses were performed in IBM SPSS Statistics 22 and R statistics. MST characteristics were compared between networks based cortical and subcortical nodes versus cortical nodes only using Wilcoxon signed ranks tests. As the correlations with multiple MST characteristics were tested, false discovery rate correction was performed to correct for increased probability of type I errors due to multiple testing.

Following the approach described by Otte et al. (2015), we used a linear mixed model to determine MST changes across adulthood. Linear or quadratic age effects were explored using fixed factors of age or age2 and sex. Akaike's information criterion was used to determine the model of best fit. For quadratic aging effects, 2000 bootstrap fits were performed to determine peak age and 95% confidence intervals.

RESULTS

Minimum spanning tree characteristics

An MSTref based on all 46 adult healthy subjects in the Netherlands dataset is shown in Figure 3. Network characteristics are presented in Table 2. We calculated which percentage of connections was the same (i.e., overlapped) for individual subjects compared to MSTref. The MST of individual subjects matched the reference MST for a mean 58.1% of connections (standard deviation (SD) 4.9%) when both cortical and subcortical regions were included, and 58.6% of connections (SD 4.7%) for cortical regions only. The network including subcortical regions had a significantly shorter diameter, lower leaf fraction, and higher kappa, tree hierarchy, and maximum betweenness centrality.

Figure 3

MST of the human brain (The Netherlands dataset). Visualization of the MST in the Netherlands dataset (group average). The left figure shows the MST for cortical and subcortical regions, the right figure shows cortical regions only [Color figure can be viewed at http://wileyonlinelibrary.com]

Table 2

MST characteristics

Measure	Reference	Individual subjects	Reference (cortical network)	Individual subjects (cortical network)	Wilcoxon signed ranks test (p value)
Utrecht dataset
Diameter	0.18	0.20 (0.03)	0.29	0.28 (0.04)	−5.90 (<.001)*
Leaf fraction	0.47	0.50 (0.04)	0.40	0.44 (0.04)	−5.65 (<.001)*
Kappa	2.74	2.92 (0.18)	2.55	2.66 (0.16)	−5.64 (<.001)*
Tree hierarchy	0.33	0.28 (0.04)	0.27	0.26 (0.04)	−3.75 (<.001)*
BCmax	0.71	0.68 (0.05)	0.75	0.66 (0.04)	−2.29 (.021)*
Degreemax	0.07	0.09 (0.02)	0.07	0.09 (0.02)	−1.27 (.206)
Australia dataset
Diameter	0.36	0.375 (0.05)	0.39	0.411 (0.05)	−11.021 (<.001)*
Leaf fraction	0.26	0.259 (0.03)	0.24	0.264 (0.03)	−4.06 (<.001)*
Kappa	2.25	2.260 (0.04)	2.23	2.264 (0.04)	−3.45 (.001)*
Tree Hierarchy	0.19	0.205 (0.03)	0.18	0.208 (0.03)	−2.70 (.005)*
BCmax	0.68	0.642 (0.04)	0.69	0.646 (0.04)	−0.94 (0.344)
Degreemax	0.04	0.045 (0.004)	0.05	0.052 (0.004)	−13.39 (<.001)*

MST characteristics for the full network and the cortical subnetwork. Measures are described for the reference network MSTref based on the mean network of all subjects, and for individual subjects mean (standard deviation). Wilcoxon signed rank tests were used to compare network characteristics of individual subjects for networks based on cortical and subcortical connections versus cortical connections only. Asterisks mark significant tests after false discovery rate correction for multiple testing.

We then tested the between‐subject variation of individual connections that are part of the MST. Figure 4 shows the MST connections weighted by the number of subjects where they were part of the MST (connections are shown when present in >25% of subjects). Figure 4 illustrates that the presence of MST connections varied across subjects. The most stable connections (i.e., connections that were part of the MST in at least 75% of subjects) are shown in Figure 5.

Figure 4

Intersubject variability of MST connections. Visualization of the occurrence of connections in the MST (MST occurrence) across subjects. Thickness of connections is higher for connections that were part of the MST in more subjects. Connections present in at least 25% of subjects are shown for clearness. The MST based on this occurrence matrix is the same as the reference MST based on the group averaged connectivity matrix, indicating that the reference is unaffected by outliers [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 5

Most stable MST connections. Twenty connections were present in at least 75% of subjects [Color figure can be viewed at http://wileyonlinelibrary.com]

The connections that were most frequently part of the MST in individual subjects were the same connections that were part of the MSTref, based on the group average connectivity matrix. To quantify this overlap, we calculated a second group average MST. In this MSToccurrence, not the group‐averaged number of streamlines was used to define connection strength, but the percentage of subjects where each connection was part of the MST of the individual subject. This MSToccurrence fully matched MSTref, that is, all connections were exactly the same.

Using FA values instead of the number of streamlines as edge weights led to a different topology of the MST (Supporting Information, Figure S1). The group average MST based on M FA showed 27.16% connection overlap with M NOS for the network including subcortical connections, and 24.69% for cortical connections only. The mean overlap of individual subjects with the reference MSTFA was 18.94% (SD 4.8%). Again, including subcortical connections resulted in a significantly lower diameter of the MST (p < 0.001; see Supporting Information, Table S1 for network characteristics of the MST based on FA connection weights). M NOS was used for further analysis, as this network was the most stable backbone across subjects.

Nodal characteristics

The degree distribution and betweenness centrality distribution across nodes is plotted in Figure 6. When subcortical regions were included in the analysis, hub nodes based on degree (mean + 1 SD) (11 in total) were the bilateral thalamus, putamen and insula, the left hippocampus, and the right caudate nucleus, superior and inferior parietal cortex, and precentral gyrus. Hub nodes based on betweenness centrality (again 11 in total) were the bilateral thalamus and precentral gyrus, left hippocampus, entorhinal cortex, parahippocampal area,temporal pole and insula, and the right putamen, and superior parietal cortex. The overlap with the rich club is visualized in Figure 7. 6/11 hub nodes in the MST overlapped with the rich club, while some others were in line with previous rich club definitions in literature (e.g., bilateral thalamus and left hippocampus). However, the definition of hubs based on the MST also showed differences with previous literature, for example, the precuneus was not identified as a hub in MST. 8/81 MST connections were rich club connections, 25/81 feeder connections, and 48/81 local connections. When the network was based on cortical regions only, the nodal degree remained relatively unchanged, but nodal betweenness centrality of most nodes increased (Figure 8), especially in the right hemisphere.

Figure 6

Node characteristics. Ranked nodal MST betweenness centrality and degree. Values represent group‐averaged means for each node [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 7

Overlap between the MST and the rich club. Rich club nodes are marked red [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 8

Effects of inclusion of subcortical regions on nodal characteristics. Effects of exclusion of subcortical regions on MST betweenness centrality and degree of cortical nodes. Color bars represent delta scores obtained by subtracting the value for the cortical MST from the value of the MST including subcortical regions. While the degree of cortical nodes remains relatively unaffected, the betweenness centrality is lower when subcortical regions are taken into account [Color figure can be viewed at http://wileyonlinelibrary.com]

ROI volume and distance effects

ROI volume was positively correlated with nodal degree (ρ = 0.48; p < 0.001) and nodal betweenness centrality (ρ = 0.42; p < 0.001). An MSTref based on number of streamlines corrected for ROI volume showed 80.5% overlap with the (ROI volume uncorrected) MSTref, indicating that 19.5% of connections are possibly influenced by ROI volume (for a visualization, see Supporting Information, Figure S2). Overlap with the reference network for the volume corrected MST was 60.2%, SD 3.6%.

Second, we calculated the Euclidian distance between the centroid voxel of each ROI as a measure of physical distance. We then plotted the occurrence of connections in the MST and the original network as a function of distance (Figure 9). The distribution of MST connections differed from the distribution of the full connectivity matrix (Kolmogorov‐Smirnov test = 0.3560; p < 0.001). The MST connections were found to be a subset of relatively short connections of the full connectivity matrix (mean Euclidean distances 25.6 and 40.8, respectively; Mann–Whitney U = −12; p < 0.001).

Figure 9

Euclidean distance plots. Euclidean distance between brain regions in general (nodal distance), the Euclidean length of tractography connections (full connectivity matrix), and the subset of connections that form the MST (minimum spanning tree). The figure illustrates that the edges that form the MST are relatively short connections [Color figure can be viewed at http://wileyonlinelibrary.com]

Replication datasets

To characterize reproducibility between different datasets and scanners, we calculated the overlap for cortical connections with a reference network based on an entirely different dataset, MSTICBM. In the Netherlands dataset, we found a mean overlap with MSTICBM for 46.5% of connections (SD 4.8%) for individual subjects, and 58.21% of connections for the group mean MST.

The reproducibility of the MST reference network was also validated in a replication dataset. The dataset contained scans of 197 subjects obtained at five scanning sites across Australia. The scans were obtained with the same scanner and analyzed with the same processing pipeline. Main differences with the Netherlands dataset were the use of a 1.5 T scanner, using different scanning protocols and preprocessing pipelines, and use of the AAL atlas instead of the Desikan‐Killiany Atlas. In the Australia dataset, MST connection overlap of individual subjects with a reference MST based on their group‐averaged connectivity matrix was 88.68%. Table 3 shows the stability of the MST for different scanning sites. The connection overlap for the group averaged MST of each scanning site was 91.01%–97.75%, suggesting that scanning site did not significantly affect MST topology.

Table 3

Dataset 3 was obtained at five scanning sites in Australia, using the same scanner and processing pipeline at each site, and was based on the AAL atlas containing 90 cortical and subcortical regions

% overlap	Average	Sydney	Melbourne	Perth	Newcastle	Brisbane
Average (N = 197)	88.68	96.63	95.51	96.63	96.63	94.68
Sydney (N = 38)	96.63	89.47	92.13	95.51	95.51	93.26
Melbourne (N = 79)	95.51	92.13	88.42	96.63	96.63	92.13
Perth (N = 29)	96.63	95.51	96.63	89.22	97.75	91.01
Newcastle (N = 18)	96.63	95.51	96.63	97.75	88.51	91.01
Brisbane (N = 33)	94.68	93.26	92.13	91.01	91.01	87.98

The table shows the percentage of MST connections that overlap between scanning sites. The MST for each site was based on the group averaged connectivity matrix of all subjects scanned on that site. The column “Average” shows the overlap with the MST based on the connectivity matrices of all subjects from the five sites. The diagonal in the table shows the mean overlap for individual subjects with the reference MST based on the group average connectivity matrix (N = 197).

In the Australia dataset, the diameter was again shorter when subcortical nodes were included in the network as compared to the cortical network only (Supporting Information, Figure S3). However, the other MST metrics showed a different pattern of effects of inclusion of subcortical nodes for the Netherlands dataset and the Australia dataset (Table 2). For example, including subcortical regions in the network led to a lower leaf fraction in the Australia dataset, but a higher leaf fraction in the Netherlands dataset.

Age and sex

Aging effects on the MST are shown in Table 4. Scatterplots showing the relation between age, age2, and MST characteristics are shown in the Appendix: results are shown for networks using the number of streamlines as edge weights and including both cortical and subcortical networks in both datasets for the whole population (Supporting Information, Figure S4), and for males and females separately (Supporting Information, Figure S5). In addition, results are shown for networks based on cortical nodes only (Supporting Information, Figures S6 and S7), and for the networks using the number of streamlines or FA as edge weights (Supporting Information, Figure S8; Netherlands data only). Kappa (a measure characterizing the homogeneity of the degree distribution in the MST) robustly correlated with ageing effects across datasets and different definitions of edges and nodes. Cortical network topology in the Australia datasets was correlated with age2 for several characteristics, while these effects were not significantly present for networks including subcortical structures. No significant gender effects were found for any of the global MST properties. Matrices representing reference MSTs for age groups 18–25, 25–35, 35–55, and 55–65 years from the Australia dataset were made available as Supporting Information. Spatial layouts of the MST of the Netherlands and Australia datasets are provided in Supporting Information, Figures S9 and S10, respectively.

Table 4

Aging effects

				Age			Age × Age
Dataset	Nodal definition	Edge definition	Metric	Regression coefficient	CI (2.5%)	CI (97.5%)	Age change	CI (2.5%)	CI (97.5%)
Australia	Cortex	nos	bcmax	‐	‐	‐	‐	‐	‐
Australia	Cortex	nos	diameter	‐	‐	‐	48.8	42.5	61.0
Australia	Cortex	nos	eccrange	‐	‐	‐	48.8	41.9	61.1
Australia	Cortex	nos	kappaa	‐	‐	‐	35.8	28.3	40.6
Australia	Cortex	nos	leafb	‐	‐	‐	35.8	2.9	39.6
Australia	Cortex	nos	tha	‐	‐	‐	38.0	32.0	42.7
Australia	Cortex and subcortex	nos	bcmax	‐	‐	‐	‐	‐	‐
Australia	Cortex and subcortex	nos	diameter	‐	‐	‐	‐	‐	‐
Australia	Cortex and subcortex	nos	eccrange	‐	‐	‐	‐	‐	‐
Australia	Cortex and subcortex	nos	kappa	‐	‐	‐	40.2	34.1	46.6
Australia	Cortex and subcortex	nos	leaf	‐	‐	‐	40.1	34.6	45.8
Australia	Cortex and subcortex	nos	th	‐	‐	‐	40.4	34.7	48.2
Netherlands	Cortex	nos	bcmax	0.067	−0.017	0.150	‐	‐	‐
Netherlands	Cortex	nos	diameter	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	nos	eccrange	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	nos	kappab	0.518	0.215	0.821	‐	‐	‐
Netherlands	Cortex	nos	leaf	‐	‐	‐	48.7	42.3	60.8
Netherlands	Cortex	nos	th	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	nos	bcmax	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	nos	diameter	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	nos	eccrange	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	nos	kappaa	0.370	0.013	0.727	‐	‐	‐
Netherlands	Cortex and subcortex	nos	leaf	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	nos	th	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	fa	bcmax	0.098	0.003	0.194	‐	‐	‐
Netherlands	Cortex	fa	diameter	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	fa	eccrange	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	fa	kappaa	‐	‐	‐	45.9	41.2	54.2
Netherlands	Cortex	fa	leaf	‐	‐	‐	‐	‐	‐
Netherlands	Cortex	fa	th	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	fa	bcmax	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	fa	diameter	−0.041	−0.092	0.010	‐	‐	‐
Netherlands	Cortex and subcortex	fa	eccrange	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	fa	kappa	0.237	−0.094	0.568	‐	‐	‐
Netherlands	Cortex and subcortex	fa	leaf	‐	‐	‐	‐	‐	‐
Netherlands	Cortex and subcortex	fa	th	‐	‐	‐	‐	‐	‐

p < 0.05.

p < 0.01.

DISCUSSION

We proposed to use MST analysis to define a reference network or standard backbone of the network of the human connectome. Using a reference network based on a group average connectivity matrix of healthy adults, we found a mean MST connection overlap of 58.1%–88.7% for individual subjects, depending on the used atlas, scanner and processing steps, and we demonstrated the reproducibility of this overlap in two large independent validation data sets. This indicates that the majority of MST connections are similar in healthy adults, but individual connections may vary between subjects. We recently showed that MST characteristics reflect properties of the underlying network (Tewarie et al., 2015), and here we tested whether the MST could be used to characterize some key anatomical aspects of the human connectome. MST hubs showed overlap with regions described in literature as hub nodes and rich club nodes in conventional graph analysis (Gong et al., 2009; van den Heuvel and Sporns, 2011).

Two recent studies showed that the MST can be used to capture biologically relevant characteristics of the human brain. Otte et al. (2015) showed that MST characteristics are sensitive to aging effects, which was replicated in the present work. A recent study in patients with psychotic symptoms showed that MST analysis is a viable method to directly compare structural brain networks of different populations, and to capture psychosis‐related alterations in brain network topology (van Dellen et al., 2015). MST analysis has also proven useful in several empirical studies on functional brain networks, showing a less integrated network in patients with neurological and psychiatric disorders including schizophrenia, epilepsy, Parkinson's disease, and multiple sclerosis (Stam, 2014).

The MST is able to characterize fundamental network properties while it remains unaffected by possible group differences in edge density or strength. In addition, the use of a reference network could also be useful for network comparison between groups, and to characterize deviations in network topology in individual subjects. Provided that similar preprocessing steps are used, the MST of individual subjects could be compared directly to the reference MST described in this study (the reference MST matrix is provided as Supporting Information).

Network specificity is another important advantage of the MST. Specificity is more desirable than sensitivity when characterizing topological properties of brain networks (Zalesky et al., 2016). The MST ensures high specificity because of the relatively sparse network representation that is achieved. For example, the MST for a 100‐node network has a connection density of 2%, whereas brain networks are typically analyzed at substantially higher connection densities (e.g., 10%–50%). Therefore, due to the higher specificity of the MST, topological properties are less likely to be influenced by the spurious connections identified with many fiber tracking methods.

The MST is acyclic by definition and may therefore not be able to fully capture network characteristics that depend on motifs and (interhemispheric) symmetry of connections. Other solutions have been proposed for the thresholding problem, such as selection of connections that are consistently present across subjects at the group level (Roberts, Perry, Roberts, Mitchell, & Breakspear, 2017). Such approaches may be appropriate when aiming to map to connectome in more detail on a group level, and take into account cyclic aspects of the network, which are by definition discarded in the MST. Rather, the MST can be used to capture a backbone of the network that reflects global characteristics of organization of the full network, including hierarchical clustering and average local clustering (Tewarie et al., 2015; Yu et al., 2015). We speculate that the MST of the connectome reflects the most important highways for information processing in the human brain. This communication likely propagates along the shortest paths between nodes or detours from these paths in the structural network (Goñi et al., 2014; Stam et al., 2015). Numerous studies on the organization of the human connectome have reported the presence of a communication backbone between cortical brain regions, consisting of a parieto‐occipital core and a high number of “rich club” connections between hubs (Gong et al., 2009; Hagmann et al., 2008; Tewarie et al., 2014; van den Heuvel et al., 2012; van den Heuvel and Sporns, 2011). We found that MST analysis is in accordance with these reports, but also that subcortical–cortical connections increase the efficiency of cortico‐cortical connectivity. As we have shown in previous work, MST hub nodes are pivotal for lowering the MST diameter and therefore network efficiency (Tewarie et al., 2015). This is in line with the description of the rich club, which also found subcortical connections to be hub nodes in our study, namely, the thalamus, putamen, and hippocampus (van den Heuvel and Sporns, 2011). Bell and Shine recently illustrated the importance of the thalamus and basal ganglia as part of a core circuit facilitating large‐scale integration of functional connectivity (Bell & Shine, 2016). In this study, subcortical connections were amongst the MST connections that were most consistently present across subjects, and these connections significantly increased the efficiency of the network. Subcortical hub regions thus play a central role in efficient communication between (cortical) brain regions, resulting in a more star‐like (i.e., more integrated) MST topology. Taken together, our findings illustrate that nodal definitions may have a large impact on the graph characteristics, and suggest that subcortical regions should be included in brain network analysis.

We aimed to further establish the biological relevance of the MST, and therefore used our data to replicate previous analyses on the relation between MST topology and age and gender. We found no gender effects on MST topology, which is in line with previous work (Otte et al., 2015). Several correlations between age and MST characteristics were found, especially kappa, and especially when cortical nodes only were considered. The global pattern that emerges from our analysis as well as previous work is that the resilience of the network to targeted attack is higher and stable during adulthood as compared to late adolescence and older age (Otte et al., 2015). While previous studies found an inverted U shape of network efficiency across the lifespan, this was not replicated with the MST approach. This finding is in line with a recent study by Otte et al. (2015), who described how cortical MST characteristics based on diffusion imaging change across the lifespan. This study together with our analysis suggests that structural brain network efficiency changes across the lifespan may be explained by global changes in connection density rather than topological reorganization. In our replication dataset of 197 subjects using the AAL atlas, we also found that including subcortical regions in the network increased network efficiency, while we found weaker correlations between MST topology and age. We conclude that not only the edge definition, but also atlas and node definition determine the sensitivity of MST characteristics to capture biologically relevant correlates. Possibly, aging affects the topology of cortico–cortical connections stronger that cortical–subcortical connections. Alternatively, the accuracy of characterization of cortical–subcortical connections may differ from cortico–cortical connections, as these connections are more centrally localized and potentially have a shorter physical distance.

Brain network analysis based on diffusion tensor imaging data comes with several methodological limitations, which may also affect the MST (Fornito et al., 2013). We found highly similar MST topologies (∼90% the same connections) for groups of subjects scanned with the same scanner at different sites. However, the use of different scanners and preprocessing methods can affect the connectivity matrix and therefore also the topology of the (reference) MST. This was shown by the considerably lower overlap of the group averaged MST based on the Netherlands dataset compared to the ICBM dataset, using the same atlas. The tractography algorithm and the definition of nodes and edge weights may all affect graph properties (Bastiani, Shah, Goebel, & Roebroeck, 2012; Fornito et al., 2013). Using FA instead of NOS for edge weights resulted in a different MST and lower overlap between subjects. The interpretation of NOS or FA in terms of connection strength is not straightforward (Jones, Knösche, & Turner, 2013); we conclude that the ranking of edge weights is more stable across subjects for NOS compared to FA. The atlas that is applied may also affect the MST connection similarity between subjects: the MST consisted of around 88% of the same connections in the Australia dataset using the AAL atlas, while this was 58.1% for the Netherlands dataset using the DK atlas. Importantly, opposite correlations between MST characteristics and age were found for the different datasets. These findings indicate that the processing pipeline and definition of nodes and edges may have a large impact on the biological interpretation of network characteristics that is not resolved by strict definitions of network topology that take the thresholding problem into account.

The MST is solely dependent on the rank of the link weights of the strongest network connections, and MST analysis may therefore be less sensitive to small differences in the signal‐to‐noise ratio between subjects. This procedure, however, may overestimate volume effects on the number of streamlines, and indeed, we found that ROI surface or volume correlated with MST characteristics. Correction for ROI volume effects affected ∼20% of connections in the MST. Secondly, long distance connections are more likely to be disturbed by random errors during fiber tracking and the number of streamlines between distant regions may therefore be underestimated. We found that the MST contains a subset of connections with a relatively short physical distance. The physical distance between ROIs was based on the Euclidean distance between the centroid voxel of each ROI, and it should be considered that the relatively short distance between subcortical regions might partially explain their stable occurrence in the MST. Steiner graphs may be used in future work to investigate the role of subcortical regions in more detail. A third potential bias when using streamline tractography is that different types of artefacts during acquisition, preprocessing and analysis may introduce distance dependent errors (Fornito et al., 2013). The effects of ROI volume, distance between nodes and basic settings for fiber tracking may interact in a complex (and unknown) manner. We therefore focused on the connectivity matrix based on the number of streamlines uncorrected for ROI volume for our main analyses. It is currently impossible to determine to what extent these correlations reflect bias of MST analysis. The MST will become a more accurate estimation of the backbone structure of the brain network with advances in fiber imaging techniques.

In conclusion, we have used MST analysis to define a reference or standard backbone network of the human connectome. We proposed a reference MST that may be used to directly quantify altered network topology in individual subjects and different populations. MST analysis is a feasible method to characterize structural brain networks that is unaffected by edge density effects. Finally, our findings provide evidence that connections between cortical and subcortical regions globally shape the backbone structure human brain, and should be included in structural brain network studies.

Additional Supporting Information may be found online in the supporting information tab for this article.

Supporting Information Figure S1

Click here for additional data file.

Supporting Information Figure S2

Click here for additional data file.

Supporting Information Figure S3

Click here for additional data file.

Supporting Information Figure S4

Click here for additional data file.

Supporting Information Figure S5

Click here for additional data file.

Supporting Information Figure S6

Click here for additional data file.

Supporting Information Figure S7

Click here for additional data file.

Supporting Information Figure S8

Click here for additional data file.

Supporting Information Figure S9

Click here for additional data file.

Supporting Information Figure S10

Click here for additional data file.

Supporting Information Table S1

Click here for additional data file.