In total, five healthy subjects participated in this pilot study (mean 37.1 ± 4.9 years std age, five males). The whole procedure involved five repeat scans for each participant 1 week apart from each other. All participants were recruited through the School of Psychology, Cardiff, Wales, UK. All participants were undergoing or had previously completed a university degree course, were right handed as assessed with the Edinburgh Handedness Inventory³ and of Caucasian origin. Exclusion criteria included a current episode or a history of neurological and psychiatric disorders, drug or alcohol abuse and medication that may have an impact on the structure of the brain. For assessment, the general health questionnaire was used (Goldberg and Huxley, 1980). All subjects provided a written informed consent.

T₁-weighted structural scans were acquired using an oblique axial, 3D fast-spoiled gradient recalled sequence (FSPGR) with the following parameters: TR = 7.9 ms, TE = 3.0 ms, inversion time = 450 ms, flip angle = 20°, 1 mm isotropic resolution, with a total acquisition time of ~7 min.

High angular resolution diffusion-weighted imaging (HARDI) data were acquired in the Cardiff University Brain Research Imaging Centre (CUBRIC) on a 3 T GE Signa HDx system (General Electric, Milwaukee, USA) using a cardiac-gated, peripherally gated twice-refocused spin-echo Echo Planar Imaging (EPI) sequence, with effective TR/TE of 15R-R intervals/87 ms. Sets of 60 contiguous 2.4 mm thick axial slices were obtained, with diffusion-sensitizing gradients applied along 30 isotropically distributed (Jones et al., 1999) gradient directions (b = 1,200 s/mm²). For further details of the MRI protocol see (Bracht et al., 2016).

Data were analyzed using Explore DTI 4.8.3 (Leemans et al., 2009). Eddy-current induced distortion and motion correction was performed using an affine registration to the non-diffusion-weighted B₀-images, with appropriate re-orienting of the encoding vectors (Leemans and Jones, 2009). Field inhomogeneities were corrected for using the approach of Wu et al. (2008). The diffusion-weighted images (DWIs) were non-linearly warped to the T₁-weighted image using the FA map, calculated from the DWIs, as a reference. Warps were computed using Elastix (Klein et al., 2010) normalized mutual information as the cost function and constraining deformations to the phase-encoding direction. The corrected DWIs were therefore transformed to the same (undistorted) space as the T₁-weighted structural images. A single diffusion tensor model was fitted to the diffusion data in order to compute quantitative parameters such as FA (Basser et al., 1994). Following the method of Pasternak et al. (2009), a correction for free water contamination of the diffusion tensor based estimates was applied (Pasternak et al., 2009; Metzler-Baddeley et al., 2012). Data quality was checked by careful visual inspection and by looking at the average residuals per DWI for each participant.

DT-MRI analysis was performed using ExploreDTI (Leemans et al., 2009) following peaks in the fiber orientation density function (fODF) reconstructed from the damped Richardson Lucy algorithm (dRL) (Dell'acqua et al., 2010; Jeurissen et al., 2013). The dRL algorithm estimates multiple fiber orientations in a single voxel and therefore provides a more accurate diffusion profile than DT-MRI-based methods estimating only one fiber orientation per voxel. For each voxel in the dataset, streamlines were initiated along any peak in the (fODF) that exceeded an amplitude of 0.05. A streamline, uniform step-size, algorithm based on that of Basser et al. (2000), but extended to multiple fiber orientations within each voxel (Jeurissen et al., 2011), was used for tractography. Each streamline continued in 0.5 mm steps following the peak in the fODF that subtended the smallest angle to the incoming trajectory. Termination criteria were an angle threshold >45° and fODF amplitude <0.05.

The automated atlas labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002) was registered to the HARDI data using a nonlinear transformation (Klein et al., 2010). The streamline termination points were coregistered to each AAL region. The numbers of streamlines connecting each pair of AAL regions were aggregated into a 90 × 90 connectivity matrix.

Connections between regions were computed by identifying the streamlines connecting each pair of gray matter ROIs. The endpoint of a streamline was considered to be the first gray matter ROI encountered when tracking from the seed location

Streamlines that did not connect to an ROI were discarded. Networks were computed for 13 different thresholds of streamline filtering by minimum contiguous length in white matter, from 0 to 6.0 mm in increments of 0.5 mm (Buchanan et al., 2014). For instance, a threshold of l mm discards any streamline that does not pass through at least l mm in white matter between gray matter ROIs.

In this section, we describe the nine adopted NWS derived from tractography.

Fractional anisotropy (FA) is calculated from the eigenvalues (λ₁, λ₂, λ₃) of the diffusion tensor. The eigenvectors (ϵ) give the orientations in which the ellipsoid has major axes and the corresponding eigenvalues give the magnitude of the peak along each axis (Basser and Pierpaoli, 1996). The mean diffusivity (MD) is the average of the three eigenvalues, while the axial and radial diffusivity are given by the largest and average of the two smallest eigenvalues, respectively (Basser et al., 1994).

The fourth NWS was based on average streamline tract length (ATL) leading to ATL-weighted networks. The fifth NWS estimated the Euclidean distance (ED) between the centroids of the two ROIs leading to the ED-weighed network. The Euclidean distance is computed in native space, so will vary across individuals.

The sixth NWS, termed streamline density (SD-weighted), records the interconnecting streamline density corrected for ROI size:

where S_ij is the set of all streamlines found between node i and node j (and S_ij = S_ji), and g_i and g_j are the number of gray matter voxels in nodes i and j. This approach leads to the construction of a SD-weighed network.

The seventh NWS is based on the volume of the tract (TV) leading to TV-weighted networks. The tract volume is computed by counting the number of voxels the streamlines of a bundle occupy and multiplying by the voxel size.

Two further NWSs were based on the number and the percentage of streamlines that connected a pair of ROIs. The number of streamlines (NSTR) is the absolute NSTR connecting two regions. The proportion of streamlines (PSTR) is the NSTR between each pair of regions, normalized to the total NSTR across the whole brain.

The adopted NWSs are called the NSTR and PSTR.

Figure 1 illustrates the nine alternative NWS and the corresponding weights from a scan of the first subject.

We integrated the different NWS via a linear integration based on the best matching of each pair of NWS-based brain networks in terms of their maximum information flow using the graph diffusion distance metric (gDDM) as described in the next section.

We topologically filtered the IWSBN using a data-driven thresholding scheme that optimizes the information flow over the cost of the surviving/selecting connections. Below, we describe the proposed data-driven topological filtering scheme.

For each of the nine weighted SBNs, we estimated six network metrics at the network level and two at the node level. Specifically, for the network level, we estimated global efficiency, local efficiency, characteristic path length, radius, eccentricity and mean weight. For the node level, we estimated global and local efficiency.

For each metric, an agreement between sessions was computed, via ICC (Shrout and Fleiss, 1979). ICC values were extracted for both network and node level and for every NWS-based brain network, for the IWSBN and also its topological filtering version IWSBN^TF.

High test-retest reliability is a prerequisite for a connectomic metric to allow for the distinguishing of different individuals (Zuo and Xing, 2014) and also for developing a biomarker of the application of functional connectomics, such as mapping growth charts of human brain function (Dosenbach et al., 2010; Castellanos et al., 2013). Therefore, beyond developing a biomarker, estimations of the test-retest reliability of functional connectomics are valuable for providing a reference regarding how strongly the estimated variables affect the observed results and guiding the significant value of the findings of both normal and abnormal brains (Zuo et al., 2014).

Recognition accuracy was assessed for each individual scan compared to the rest based on a k-nearest neighbor (k-NN) classifier with k = 4 and adopting a leave-one-out cross-validation scheme (LOOCV). Instead of the Euclidean Distance (ED) most commonly used in k-NN classifiers, here we used the proposed gDDM (see section Graph Diffusion Distance Metric). gDDM is a more appropriate metric compared to ED to quantify the distance between two SBNs regarding their distance in terms of information flow based on their topology. The proposed gDDM metric is based on the eigenanalysis of the Laplacian matrices with known attributes in terms of graph theory and diffusion processes (Fouss et al., 2012).

To demonstrate the effectiveness of the proposed method in a binary classification problem, we analyzed a large dataset consisting of 123 individuals with psychotic experience and 125 age and gender-matched controls. The details of the cohort and the MRI scanning protocol can be found in the original publication (Drakesmith et al., 2015a).

We followed a binary classification procedure with a 10-fold cross-validation, employing as an input, each weighted SBN separately but also the IWSBN and the topologically filtered version (IWSBN^TF). As a classifier, we used a tensor subspace analysis to reduce the initial high-dimensionality of the original functional connectivity network to a space of condensed descriptive power (Dimitriadis et al., 2015b,c; Antonakakis et al., 2016). The input on TSA is a 3D tensor-matrix of dimensions (subjects × ROIs × ROIs). As a classifier, we used a support vector machine with RBF kernel.

The proposed IWSBN^TF was derived after first linearly combining the nine NWS and after that, topologically filtering the outcome IWSBN. Our second thought was to weight each node independently within each of the nine NWS first and secondly to weight the whole NWS-network with the proposed methodology. To get the linear weights l_w for each node within each NWS-network, we followed five strategic network lesion schemes, three based node-wise and two cluster-wise.

The first three node-wise strategies were: (1) zeroing half of the connections of each node; (2) diminishing the weights of the connections of each node by 50%; and (3) combining both where half of the connections of each node were zeroed while the weights of the second half were diminished by 50%. The three node-wise attack strategies were followed for one by one node at every NWS and then we estimated the gDDM distance between the original network and the attacked network. Finally, the derived vector with the 90 gDDM values was normalized such as its sum was equal to one. Then, we multiplied each NWS-network node-wise with this vector and afterward with the network-wise approaches based on the present methodology.

The two cluster-wise lesions were: (1) the distinction of the whole set of hubs into rich-club hubs and non-rich-club hubs (van den Heuvel and Sporns, 2011); and (2) the functional clustering of the NWS-network into distinct clusters using the modularity algorithm (Newman, 2006). The two cluster-wise attack strategies were followed for each cluster at every NWS, based on the three node-wise attack strategies targeting either connections within rc-hub subgraphs and/or non-rc-hub subgraph connections and/or the connections between the non-rc and rc-hubs. Then, we estimated the gDDM distance between the original network and the attacked network. Finally, the derived vector with n_cluster gDDM values was normalized such as its sum was equal to one.

The whole procedure was added as a first step before the proposed network-wise linear combination of the NWS-network into a single IWSBN All the NWS-networks were pre-filtered with the proposed data-drive thresholding scheme. The node-wise linear weighting step prior to the proposed network-wise was evaluated based on the ICC values of the adopted network metrics both network and node-wise. Additionally, the recognition accuracy of each subject scan over the rest of cohort was compared to the proposed method.

Dissimilarity matrices (DM) based on each NWS across scans and subjects were estimated based on the gDDM. Figure 5 demonstrates the DM for each of the NWS across repeat scans and the related graph embedding based on multidimensional scaling (MDS). The NSTR proved to better discriminate the five subjects compared to the rest of the methods.

The ICC scores were excellent—ranging from 0.75 to 1—for six out of nine network metrics for the entire set of network metrics. These NWS include the ATL, BIN, SD, ED, NSTR, PSTR, and TV (Figure 6). The related group-averaged values of the network metrics for each NWS are shown in Figure 7. The proposed IWSBN yielded good ICC values but these were lower than those obtained for each of the six NWS (Figure 8A). Significantly, we observed an improvement of the ICC on the IWSBN^TF which reached the level of the six best NWS in terms of ICC scoring (Figure 8B). Figure 9 demonstrates the group-averaged values of network metrics on the network level.

The analysis of ICCs on global and local efficiency node-wise on the nine NWS and in both IWSBN and IWSBN^TF revealed important trends. Firstly, the ICC values for each of the NWS failed to reach a fair value (ICC < 0.1). Secondly, the ICC values derived from the IWSBN showed a large variability but reached on average ICC = 0.68 ± 0.10 for global efficiency and ICC = 0.68 ± 0.17 for local efficiency (Figure 10A). Third, the ICC values for both network metrics were improved in IWSBN^TF, reaching a mean ICC = 0.75 ± 0.02 for global efficiency and a mean ICC = 0.84 ± 0.02 for local efficiency (Figure 10B). Applying a Wilcoxon Rank Sum Test between the two distributions for each network metrics, we observed significant improvement of ICC values for IWSBN^TF (global efficiency: p = 0.0035 × 10⁻⁷, local efficiency: p = 0.0067 × 10⁻¹⁰). Figure 11 demonstrates the group-averaged values of network metrics on the node level.

Dissimilarity matrices (DM) based on both IWSBN and IWSBN^TF across scans and subjects were estimated based on the gDDM. Figure 12 demonstrates the DM for both IWSBN and IWSBN^TF across repeat scans and the related graph embedding based on multidimensional scaling (MDS). The proposed topological filtering scheme improved the discrimination of the five subjects compared to the original IWSBN.

Applying a k-NN classifier with k = 4 and gDDM as the appropriate distance metric under a LOOCV scheme, we succeeded to accurately classify each scan to the right person based on IWSBN^TF. Similar results were also obtained with NSTR (Figure 5).

For a better estimation of the discrimination of proposed IWSBN^TF with NSTR, we defined the following quality of clustering index (QCI) (Dimitriadis et al., 2012):

The QCI quantifies the average similarity of the IWSBN^TF across scans within each subject (cluster) expressed in the denominator of Equation (6) and the average dissimilarity between every pair of subjects (clusters) across their scans expressed in the numerator. Both similarity and dissimilarity were estimated based on the gDDM. The higher the dissimilarity between the clusters (numerator) and/or the lower the dissimilarity within the clusters (denominator), the higher the QCI. The numerator is averaged across all possible combinations of subjects (clusters) while the denominator across subjects (clusters). The first term was used to equalize the effect of between-subject (clusters) comparisons vs. within-subject comparisons (clusters). This inversed coefficient guarantees that in the case of all the weights in the DM being equal then it will take a value of one. Therefore, the higher the value of the QCI above one, the higher the separability between the network topologies of the subjects.

We first tabulated all the IWSBN^TF across scans and subjects into a 4D graph with dimensions equal to [subjects × scans × Rois × Rois] called D_ IWSBN^TF. Afterward, we estimated the QCI for each NWS and for both IWSBN and IWSBN^TF.

The QCI was 1.45 for IWSBN and 6.45 for IWSBN^TF while for NSTR the QCI was 4.57. This result can be interpreted as a higher separation of network topologies with our approach compared to the best of NWS.

To further enhance the integration of the nine alternative WNS for the construction of an integrated SBN, we repeated the whole procedure splitting the nine weighted versions of SBN into three triads (the first three, the second three and the last three). Figures 13–15 illustrate the DM and their embedding into a 3D-space. The highest separability between the network topologies of the subjects have been demonstrated for ATL, SD, and ED, while the worst for NSTR, PSTR, and TV, where three subjects overlapped on the same embedding space (Figure 15B). The QCI score was lower compared to the original approach where we combined the nine weighted SBN (Figure 12).

Our classification results demonstrated a higher classification accuracy (65.3%) between the two groups for the proposed IWSBN^TF. The classification performance of the nine weighted strategies was lower than by chance (<50%; see Table 1). Additionally, the data-driven topological filtering via the OMST algorithm (Dimitriadis et al., 2017a,b) further improved the classification accuracy (from 57.23 to 65.34%; see Table 1).

Analysis of reliability is challenging for neuroimaging as the main results presented in a study depend on both the adopted network metrics and the metrics used to characterize the weight of a connection. Additionally, many neuroimaging studies based on SBNs attempt to shed light on developmental differences, differences between clinical populations and also between a control group and a disease group. A basic reason why all these proposed connectomic biomarkers are not used in daily practice in hospitals is their reliability (Dimitriadis et al., 2015a). A second reason is that there are many studies on the same topic (e.g., brain disease) based on small datasets that adopted different NWS and arbitrary topological filtering schemes. Meaningful aggregation of these, even in the case that their metadata are free available, is impossible. A third reason is that until now, it is not standard practice to assess the reliability of network metrics derived from SBNs across repeat scans on the same population but with different dMRI protocols and scanners (3T vs. 3T or 3T vs. 7T). This is an issue that we would like to investigate in future studies with the proposed scheme.

A basic issue with the ICC is that it needs large samples in order to estimate scores to acceptable precision. A study estimated that for two repeated measures, in order get an acceptable ICC score of 0.8 with a 95% confidence interval of 0.2 width, then at least 52 subjects are needed (see Table 3 in Shoukri et al., 2004). In a similar vein, for an ICC score of 0.6 with 95% confidence intervals of 0.2 width, repeated measures from 158 subjects would be required. Clearly for most MRI-based studies, this scenario of repeated scans for hundreds of scans is unrealistic. However, our analysis provided a methodology of how to combine different NWSs into a single integrated graph based on a gDDM that counts the network topology as a whole and quantifies the distance between two SBNs in terms of their information flow. The results of the proposed methodology presented here, even in a small dataset, are of paramount importance since they are completely data-driven in any preprocessing step. Simultaneously, they provided novel directions of how to untangle hidden information within dMRI-based brain networks by working on an individualized manner and without any averaging approach (group or scan-wise). Additionally, we proposed a data-driven thresholding scheme applied to the IWSBN that improved the ICCs of basic network metrics at both node and network levels compared to the original IWSBN and each of the adopted NWS. Our data-driven method can be seen as a methodology for improving network reliability on SBNs (Zalesky et al., 2010; Drakesmith et al., 2015b). Complementarily, our approach provided excellent discrimination of the network patterns of the five subjects based on the recognition of each scan to the targeted subject. Finally, our method better separated the five groups of scans based on the topological filtering version of IWSBN compared to the best NWS.

It is important to mention here the limitations of the current study due to the small dataset for exploring the test-retest reliability statistics. In the era of open science resource where multisite worldwide neuroimaging labs share neuroimaging datasets, it is important to demonstrate novel techniques that improve the reliability of connectomics in common neuroimaging data (Zuo et al., 2014). A recent Consortium for Reliability and Reproducibility (CoRR) is working to address this gap and establish test-retest reliability as a minimum standard for methods development in functional connectomics (Zuo and Xing, 2014; Zuo et al., 2014) and morphological measurements (MacLaren et al., 2014). Reliability is important to build reliable connectomic biomarkers across multi-site (Nielsen et al., 2013; Abraham et al., 2017) and also longitudinal trajectories of structural and functional brain networks across the life-span (Zuo et al., 2017).

Our future goal is to test the proposed methodology in a larger sample for validating the test-retest reliability of our scheme and also on multi-site diffusion-based structural brain networks for building reliable connectomic biomarkers.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.