NBS-SNI, an extension of the network-based statistic: Abnormal functional connections between important structural actors

Shedding light on the heterogeneous structure-function relationships in the brain has been the aim of many studies over the past 15 years (Baum et al., 2017, 2020; Honey, Thivierge, & Sporns, 2010; Park & Friston, 2013; Pessoa, 2014; Suárez, Markello, Betzel, & Misic, 2020; Witvliet et al., 2021). These endeavors stem from the well-established hypothesis in neuroscience that functional activity in the brain should be somehow related to its underlying structural architecture; indeed, the structural connectome is what supports the vast collection of functional patterns available to the human brain (Sporns & Kötter, 2004).

Correlations between brain structure and brain function have previously been demonstrated to be altered in schizophrenia, autism, bipolar disorder, small vessel disease, and mild cognitive impairment (Cao et al., 2020; Cocchi et al., 2014; Rudie et al., 2013; Tay et al., 2023; van den Heuvel et al., 2013; Wang et al., 2018; Zhang et al., 2019). This structure-function coupling has also been shown to be prominent during youth, which would make it a crucial period for proper cognitive maturation (Baum et al., 2020).

However, integrated frameworks for concurrently analyzing both structure and function remain limited. Prior work has used disparate methodological approaches for relating brain structure and function, including the rich-club coefficient (van den Heuvel et al., 2013), individual or group level SC-FC correlational analyses (Cao et al., 2020; Tay et al., 2023; Wang et al., 2018), or primary component analysis (Rudie et al., 2013). These approaches may disregard information shared between structural and functional representations as structural and functional properties are assessed separately, and then combined. This may also result in confounds when comparing studies, due to the different statistical measure or metrics used.

One popular tool for evaluating differences in edge properties between (case and control) groups is the network-based statistic (Zalesky, Fornito, & Bullmore, 2010). This approach has previously shown group level differences in functional connectivity edge weights in disorders including schizophrenia (Zalesky et al., 2010), autism (Pascual-Belda, Díaz-Parra, & Moratal, 2018), bipolar disorder (Roberts et al., 2017), and Alzheimer’s disease (Zhan et al., 2016). However, one key weakness of the NBS is that it only assesses properties of edges, and only assesses one edge property at a time. We look to extend the previous work undertaken in this direction by creating a framework that combines information from multiple streams, incorporating differences in both nodes and edges between groups.

Here, we introduce a procedure called Network-based statistic–simultaneous node investigation (NBS-SNI), which assesses properties of both nodes and edges, allowing for simultaneous assessment of two representations of the brain (e.g., structural and functional networks). In typical case-control paradigms, NBS-SNI aims to identify differences between two groups (case/control) of brain networks through an examination of abnormal connectivity between the most important actors of these networks. It does so by integrating two network representations of the brain at once, in pursuance of a more comprehensive representation of the brain. Within the proposed NBS-SNI framework, two properties are computed in each network of the healthy and the condition groups (one for each edge and one for each node), and compared through the lens of statistical testing. Combining information from nodes and edges in this way may increase both pathophysiological understanding and statistical power.

In the Methods section, we give an overview of the general framework of the method and discuss potential choices of node properties that could be employed with our approach. In the section Contrast Identification, we test our method with synthetic data to identify the regimes where our method NBS-SNI offers a statistical advantage over NBS or a generic link-based method. By considering the principal actors of the structural networks, we show that, in certain regimes, our method can flush out false positives and increase the resolution at which significant edges (true positives) can be identified. Finally, in the section Brain Network Classification Performances, we apply our method on two real case-control datasets and perform classification of brain networks using various node properties. As an illustration, in the first dataset, using brain networks derived from neuroimaging data of an autism-healthy study (Rudie et al., 2013), where both the structural and functional networks are available (N = 29 ASD, N = 19 typically developing; TD), we found that NBS-SNI offers an improvement of 9 percentage points in prediction performance (73%) compared to the standard NBS (64%). On the second dataset, comprised of early psychosis-healthy individuals (N = 75 EP, N = 45 typically developing; TD), we incorporated anatomically derived measures (intrinsic curvature index and gray matter volume) as node properties in our proposed framework. We found that NBS-SNI offers a 4.9 percentage points gain compared to NBS.

NBS-SNI is an extension to the well-established network-based statistic (NBS; Zalesky et al., 2010), a general framework used to test specific hypotheses in typical case-control studies. In essence, the formulation of NBS-SNI we propose tries to identify the most important nodes, common to the healthy and condition brain networks, and probes for abnormalities (hyper or hypo) in the interconnection around these influential brain regions. This is different from the standard NBS, which considers the extent to which significant edges (from a t test perspective) are interconnected to form a component. The notion of component, that is, retaining only the significant edges that form a connected component, is relaxed with this approach, since an additional restriction comes from the consideration of significant nodes, again by means of a statistical test, in addition to significant edges. A node property is computed for every node in the brain networks. The choice of a relevant node property in the context of brain network analysis is at the liberty of the investigator; a few possibly relevant node measures for network neuroscience are presented in the Methods section.

It is hypothesized that such a set of edges where stronger connections are observed between particular nodes in one group versus the other could be associated with certain symptoms, or network pathologies. Therefore, the choice of an appropriate node measure can be guided by prior knowledge of the specific condition under investigation. The resulting set of connections identified by this approach can then be paired with the framework NBS-predict proposed by Serin, Zalesky, Matory, Walter, and Kruschwitz (2021) to yield a prediction of individuals’ classes in case-control studies. In this article, the relevance (significance) of the subnetworks extracted from real data are assessed by the prediction performances that they yielded on individual brain networks through the use of NBS-predict framework (as opposed to significance being assessed by permutation testing, as it was originally done in Zalesky et al., 2010).

This approach allows the analysis of both structural and functional representations of individual brain networks at once, under a single framework. Owing to the statistical framework of the method, which compares individuals across two representations of brain networks, a perfect node correspondence is required across both representations.

Depending on the imaging modalities employed, the nodes representing the constituents of the brain network representation could be individual neurons (microscale), groups of neurons (mesoscale), or brain regions (large-scale). The interactions (or associations) between the nodes are represented by the edges of the network, and their nature depends on the imaging technique employed. A brain network of n nodes is represented by a n × n connectivity (or adjacency) matrix, where the entries represent the pairwise associations between every possible pair of nodes.

The real brain networks investigated in this work originate from both fMRI (functional) and DTI (structural) recordings (large-scale) for every individual; micro/mesoscale brain networks will not be discussed here. Interactions measured by DTI structural imaging are bidirectional (undirected) edges whose weight corresponds to the density of white matter fibers connecting two brain regions. Interactions measured by functional imaging are also weighted and bidirectional, but their weight represents a measure of temporal correlation between the recorded dynamic activities between each pairs of parcellated brain regions.

NBS-SNI requires the use of a node property, which should be chosen to capture relevant node features related to efficient communication, resilience, clustering, or to the symptoms associated with the specific condition under investigation. Some of these node properties could be the degree, or the betweenness/closeness centrality, as they are well suited for a large-scale brain networks dataset, in which the correspondence between the parcellated brain regions (the nodes) is reliable and the volume of each parcellated region does not vary much across individuals (Simpson, Lyday, Hayasaka, Marsh, & Laurienti, 2013). Node features could also be anatomically derived measures, such as the cortical thickness, gray matter volume, curvature, etc. The chosen node property is computed for every node in each network of the healthy and the condition groups, and then compared through the lens of statistical testing. Owing to the central limit theorem, it is worth mentioning that due to the smaller number of node properties (n) compared to the number of edge weights (n(n − 1)/2), the normality assumption may not always be reasonable for the computed node properties. In such a case, nonparametric statistical tests might be considered, instead of the usual Student’s t test, which will be used for edge weights statistical testing. Before getting into the details of these statistical tests, we briefly describe the node properties used here, and explain their relation to the notions of efficiency and resilience in brain networks. We also want to highlight the fact that the node properties don’t necessarily have to be computed from the structural brain networks. They could also be computed directly from the functional networks themselves. Moreover, in the second real dataset that we test with our toolbox, we used anatomically derived node measures from T1 images as a node property.

The interpretation of closeness centrality is usually one of access of efficiency or independence from potential control by intermediaries (Brandes, Borgatti, & Freeman, 2016). From this interpretation, a node of high closeness centrality would be less vulnerable to potential failures of other brain regions to which it is connected, thereby easing efficient rerouting of information along different pathways. In the context of case-control studies of brain networks, this interpretation of the closeness centrality makes it an interesting property to identify important nodes or identify between-group abnormalities. Closeness centrality will be used for both synthetic (see section Contrast Identification) and real data (see section Brain Network Classification Performances).

We also consider another measure called the current flow closeness centrality (sometimes called information centrality; Stephenson & Zelen, 1989). Contrary to the closeness, this node property does not make the assumption that communication happens solely along shortest paths (geodesics). Instead, information centrality uses all paths, and assigns them a weight based on the information they contain, therefore potentially capturing more subtle effects caused by nongeodesic paths in complex systems. This distinction with the closeness centrality makes it an interesting nodal property to investigate with NBS-SNI’s framework. The reader is referred to Stephenson and Zelen (1989) for detailed calculations of this measure. This measure will be used in section Brain Network Classification Performances.

The last node property that will be investigated with NBS-SNI in section Brain Network Classification Performances with real data is called the local information dimension (LID) (Wen & Deng, 2020), a property derived from the local dimension (LD) (Silva & da F. Costa, 2013). For brain networks, a local definition of dimensionality allows one to quantify the local resilience of the network. The logic behind this statement is that objects of higher dimensionality might be less susceptible to being compromised by a random attack, since higher dimensionality often entails more redundant pathways, that is, alternative information transmission pathways.

The LID D_{I_v} of a node v is obtained from the slope fitting of I_v(l) versus ln l. In contrast to the LD, the size of the box l with LID varies from one to half of the maximum shortest distance from the central node v. This reduces the computational complexity and focuses on the quasilocal organization around the central node v. The LID was chosen over the LD as a node property to investigate in this article owing to its reduced computational complexity, its greater effectiveness, and reasonability at identifying influencers in complex networks (Wen & Deng, 2020). Additionally, we found that using the LID as a node property resulted in more accurate predictions than the (more global) local dimension (LD) (Silva & da F. Costa, 2013).

In section Brain Network Classification Performances, we extend the LID to weighted networks by discretizing the l space into an arbitrary number of steps, starting from the minimum weighted shortest distance to half the maximum weighted shortest distance from the central node.

Finally, the gray matter volume of each parcellated brain region was also employed as a node property with the HCP early psychosis dataset in section Brain Network Classification Performances. Gray matter volume reductions and abnormalities have been reported in patients with schizophrenia (Gur et al., 2000; Gur, Turetsky, Bilker, & Gur, 1999; Hulshoff Pol et al., 2002; Vita, De Peri, Deste, & Sacchetti, 2012) and also in patients with autism (Dekhil et al., 2019), thereby making it a suitable node property to incorporate into NBS-SNI.

In the following description, the method is designed to identify the most important nodes (brain regions) common to both the condition and the control group in the structural representation, and probe for hyperconnected (or hypoconnected) functional connections between these brain regions.

In what follows, both the structural and functional connectivity matrices are available for each individual in both groups (case and control). Moreover, the method assumes that there is a “perfect” correspondence between the nodes across the networks for both representations (node-aligned networks). Nevertheless, we want to point out that it would also be reasonable to compute the node properties directly from the functional connectivity matrices. Furthermore, node properties could also be derived from anatomical measurements, without necessarily having to be extracted from the structural connectivity matrices. This will be illustrated in section Brain Network Classification Performances, where the anatomically derived measures of intrinsic curvature index and gray matter volume will be used as a node property with NBS-SNI.

Let there be N₁ pairs of networks (functional and structural) of n nodes in group 1 (condition) and N₂ pairs of networks of n nodes in group 2 (control). To find significant nodes common to both the condition and control group, a node property is computed for each node from the structural connectivity matrices. This results in a N₁ × n matrix where each row contains the node measures of all the nodes in a structural network from group 1. The same procedure is done with the other group to obtain a N₂ × n matrix of node properties (step 1 on Figure 2). Two independent one-sample statistical tests are then performed on this property, yielding two vectors of length n containing the statistic of each node, that is, one vector for each group (step 2). It is important to note at this point that the normality assumption required to apply the t test might not always be met when using node properties due to the small number of nodes (compared to the number of edges). In that case, the user might resort to a nonparametric test such as the Kolmogorov-Smirnov test instead of the one-sample t test.

A threshold t_node is then applied to the two resulting statistics independently in each group to identify significant nodes, and the intersection between the significant nodes of both groups is taken (step 3). These suprathreshold nodes correspond to the most important actors (nodes) of these networks mentioned earlier. Alternatively, it would be reasonable to retain only the suprathreshold nodes from one of the two groups if the intersection between the two sets is empty, or if one set contains many suprathreshold nodes, while the other contains a few.

In parallel, a two-sample t test is also performed on the edge weights of the functional connectivity matrices for a between-group difference (step 4), resulting in a single n × n matrix that contains the t statistic of each edge weight as its entries. A threshold is then applied to this matrix (step 5), resulting in a set of suprathreshold edges (step 6).

Finally, the intersection between the resulting suprathreshold nodes and suprathreshold edges is taken (steps 7 and 8). Three cases naturally arise. The most restrictive case consist in retaining only the edges that are shared between the “important” nodes (d = 0). Note that for d = 0, the subnetwork yielded by NBS-SNI is always a connected component. The second case implies retaining all edges that connect at least one suprathreshold node (d = 1). The third case retains edges that connect nodes that are at most at a distance d > 1 from the suprathreshold nodes. Note that this third case requires computing shortest distances from the central nodes, which increases the computational complexity of the method. Hereinafter, we chose to perform the statistical test on the edge weights of the functional connectivity matrices, while the statistical test on the node properties takes place on the structural connectivity matrices, but swapping these choices could be valid as well. Also, note that an alternative formulation where the investigator first probes for abnormalities in the nodes from the structural representation could also be valid (i.e., performing a two-sample statistical test for between-group differences), as opposed to finding nodes considered important and common to both the condition and control group. This formulation employing a two-sample t test will be used with a real dataset in section Brain Network Classification Performances. using the anatomical measure of gray matter volume as a node property on the HCP early psychosis dataset.

In our simulations presented in Methods section, NBS and NBS-SNI are sometimes compared against the false discovery rate (FDR), that is, the expected ratio of false positives among all the positives, which serves as a link-based controlling procedure. With FDR, every edge’s t score is treated independently, as opposed to NBS or NBS-SNI, which do not yield a collection of individual significant edges, but rather a connected component that they form.

We adapted the framework NBS-predict (Serin et al., 2021) to investigate the impact information about nodes taken into account by NBS-SNI has on the accuracy of a classification task. In this context, individual brain networks (functional connectivity matrices) are classified (diagnosed) into either the control or the condition group. The edge weights of the subnetworks extracted by the method are the features that are fed to various machine learning (ML) algorithms (the features associated to the edges that were not selected are set to 0). A K-fold cross validation procedure is repeated r times (see Figure 3). Finally, the method outputs an average prediction accuracy for the brain networks and a weighted adjacency matrix, in which the weights of the edges represent their individual prediction accuracy across folds. Figure 3 illustrates the workflow of the classification procedure NBS-predict.

To visualize the contribution of each edge to the overall prediction performance, the selected edges in the extracted subnetworks across folds are constantly updated and assigned a weight (encoded in a single weighted adjacency matrix), which corresponds to the prediction performance at that specific fold (here again, edges that are not selected are assigned a weight of 0). In the end, all the weights are averaged and rescaled so that the maximum prediction weight is 1. One may then visualize the selected edges across folds that contributed the most to the prediction (i.e., the largest weights denote the edges that had a greater prediction accuracy and that were selected the most across folds).

The classification ML algorithms employed in this work are Ridge Classifier, Random Forest Classifier, Decision Tree, Support Vector Machine (SVM) and Logistic Regression, all implemented from the scikit-learn (Pedregosa et al., 2011) Python package.

To compare the statistical resolution of NBS-SNI with NBS to identify a contrast, synthetic data was created for both structural and functional connectivity brain networks in a typical case-control study setup. As mentioned previously, we expect NBS-SNI to draw its strength from considering both representations simultaneously, when an interplay of contrasts exists across the two representations.

By increasing the connection probability between u and v if either belongs to n_c, the last term in Equation 5 effectively increases the closeness centrality of the nodes in n_c compared to the other nodes. In accordance with the formulation of NBS-SNI we propose in this article, both groups of synthetic brain networks (condition and control) were induced the same high centrality nodes, that is, the high centrality nodes in n_c were common to both groups of networks. The nodes in n_c serve as the important actors previously mentioned. Note that all connections are bidirectional, that is, when a connection is added from node u to node v, the same connection is also added from node v to u. Note also, that an edge can only be added once, therefore the probabilities are adjusted accordingly during the generative process.

To obtain the parameters associated with E(u, v), we plotted the probability of connection P_uv between node u and v for the real structural network dataset of 264 brain regions (UCLA Autism; Rudie et al., 2013) as a function of the distance separating them d_uv, considering both groups (ASD and control networks) altogether. We found that this relationship could be closely approximated by three different functions according to the range of distances. For distances ranging from d_uv = 0 up to d_uv = 18 mm, a function a of the form P_uv ∝ A₁$logduv−bduv−b$ was fitted. The second part was approximated using a function of the form P_uv ∝ A₂e^−ℓd_uv(d_uv)^−η₂ for distances ranging from 18 ≤ d_uv < 39 mm. For the last part (d_uv ≥ 39), a power law of the form P_uv ∝ A₃(d_uv)^η₃ was fitted.

Since structural brain networks are generally represented as weighted networks, the distribution of weights from the real autism structural dataset was also extracted. We found that a function of the form p_w ≈ α_w ⋅ e^(kw)/w + C_w, where p_w is the probability of finding a connection with weight w between two brain regions was a good approximation of the first half of this relationship. The second half (the right tail) of this relationship was approximated by a power law distribution ∝ A_w ⋅ w^−η_w, which allows for a few larger connection weights that cannot be accounted for by the first function (A_w ⋅ e^(kw)/w).

For E(u, v), the parameters we found empirically are (first part; A₁ = 133.49 and b = −899.3), (second part; A₂ = 0.000943, ℓ = 0.2277 and η₂ = −3.83), and (third part; A₃ = 275.178 and η₃ = 2.031). For the weights, we found α_w = 0.194, k = 1.038 × 10⁻⁶, and C_w = −0.0007 for the first part, and A_w = 627.78 and η_w = 2.99 for the second part. Moreover, standard deviation σ_nc = 2 was kept constant for all our experiments.

The synthetic functional brain networks are fully connected undirected weighted networks, that is, a bidirectional weighted connection exists between every pair of nodes (u, v).

In the other group (noise only group) all the entries in the adjacency matrices are sampled from 𝒩(0, σ). For σ = 1 (the value employed for all the functional contrasts in our simulations), this procedure has the effect of creating a contrast-to-noise ratio CNR = $μei$ on the edges in $μei$ in one group, while the other group consists of adjacency matrices with all their entries sampled from Gaussian noise of zero mean and variance σ².

The synthetic networks for this experiment were built using the same nodes (n = 264; including their positions) as in the UCLA Autism dataset (Rudie et al., 2013), of which ∣n_c∣ = 5 were chosen as being more central (see Figure 4A). The same number of weighted edges, 3,188, were randomly added using Equation 5, with parameters μ_nc = 1.5 and σ_nc = 2. We generated 60 synthetic structural networks equally divided into the condition and the control groups, each network consisting in a single connected component.

To build the functional networks from the condition group, a first set of eight edges connected to at least one node in n_c were given a weight using Equation 7 with $μe1$ = 0.6 and σ = 1 (subcontrast edges; dashed lines in Figure 4B). A second set of 11 edges located in the periphery of the 5 nodes in n_c were given a weight using Equation 7 with $μe2$ = 0.8 and σ = 1 (contrast edges; solid lines in Figure 4B). All other edges, as well as all edges in the functional networks belonging to the control group, were given a weight according to Equation 7 (sampled from 𝒩(0, 1)).

Our objective here is to illustrate how taking into account information about the nodes (here present in the structural networks) can allow the identification of edges that do not necessarily have the largest contrast-to-noise ratio (these subcontrast edges are shown in Figures 4C and 5A). Figure 5 shows a comparison of the performances between NBS and NBS-SNI with d = 1.

Blue edges correspond to connections identified by either method that did not belong to the subcontrast (false positives), while red edges correspond to connections that were missed (false negatives). Figure 5B shows the results for two sets of parameter thresholds t_edge obtained with NBS-SNI, which perfectly recovers the subcontrast edges. Figure 5C shows that t_edge could not be tuned for NBS to perfectly recover the subcontrast edges. For NBS, the threshold t_edge was increased to a larger value of 4.7 compared to NBS-SNI, where t_edge = 4. This is due to the fact that NBS-SNI is more restrictive than NBS, owing to the statistical test on the nodes. For NBS, this figure shows that past a threshold of 4.7, the method discards edges that belong to the contrast, even though it still misses some true positives.

Figure 5 illustrates what was already expected owing to the design of NBS-SNI, that is, it can recover a subcontrast that is embedded inside another contrast of greater CNR, a task NBS is not suited for. Evidently, this conclusion is only valid if the nodes in the subcontrast share some common property (e.g., the closeness centrality in this example), which distinguishes them from the other nodes. In our example, this shared property between the nodes in the subcontrast was induced in the structural network models, but a nodal property could also be computed from the functional network representations. We chose this formulation since nodal properties based on shortest distances seem most suited for structural networks than for functional ones.

The next experiment aims at getting a better grasp of the cases when NBS-SNI offers a gain over NBS (and when it does not). Alongside NBS-SNI and NBS, we also consider the false discovery rate (FDR) to compare the two methods with a link-based controlling procedure.

Indeed, FDR treats every edge independently, as opposed to NBS and NBS-SNI, which cannot declare individual edges to be significant, but rather the component that they form. We evaluated the specificity and sensitivity of both NBS-SNI and NBS, along with the false discovery rate (FDR).

The experiment goes as follows. This procedure is repeated 1,000 times for each set of parameters and thresholds (see Figure 6 for details). Again, synthetic networks of 3,188 edges were generated using the 264 nodes and their pairs of distances from the UCLA Autism dataset (Rudie et al., 2013). Among these nodes, ∣n_c∣ = 5 are randomly selected to be the central nodes in the structural network. Among the $∣nc∣2$ edges between these nodes and two other nodes (not present in n_c) in the functional networks, a component of 20 edges consisting of a single contrast with mean CNR $μe1$ and variance σ = 1 is selected as the set $ne1$⁠. The component of functional edges therefore encompasses the five central nodes, plus two other nodes which do not possess higher induced centrality. Using Equations 5 and 7, 60 synthetic networks equally divided into the condition and the control groups are generated. (Recall that edges in the set $ne1$ in the functional networks associated with the control group are no different than the other edges.) The most significant edges are then identified using the three methods (using d = 1 for NBS-SNI), and the true positive (TPR) and false positive rates (FPR) are computed. For NBS-SNI, d = 1 was used since some of the 20 connections in the contrast reach two nodes that do not belong to n_c.

The results are summarized by means of the receiver operating characteristic (ROC) curves shown in Figure 6. In all cases, the value of the parameter thresholds ranged from t_edge = 0.5 to t_edge = 6 (for NBS-SNI and NBS) and q = 0.001 to q = 0.99 (for FDR; q-values represent the expected false discovery rate), with 100 equally spaced values in between. For NBS-SNI, the parameter threshold for the nodes was kept constant at t_node = 0.5 for all realizations. For the synthetic functional networks, all the edge weights were sampled from a Gaussian distribution with with μ = 0 and σ = 1.

In Figure 6A–C, the parameter μ_nc in K_nc(u, v) ∼ 𝒩(μ_nc, σ_nc = 2) was varied in the structural network models for a given CNR = 0.3 (⁠$μe1$ = 0.3; σ = 1) between the two groups in the functional network models. Since only the properties of the structural models are altered, NBS and FDR’s performances remain the same across Figure 6A–C. It is only the performance of NBS-SNI that improves, that is, the curves’ ascension is steeper, as the parameter μ_nc in the structural model is increased.

Figure 6D–F shows the ROC curves for three different CNRs in the functional network models, for a constant parameter μ_nc = 1.5 in the structural network models. The gain in statistical resolution offered by NBS-SNI compared to the others becomes less significant as the CNR is increased.

In Figure 6G–I, the CNR was again varied, but for a constant structural model parameter increased to μ_nc = 2. On the right, Figure 6F with CNR = 1, all three methods offer similar performances after FPR = 0.01, thereby showing the gain in statistical resolution of NBS-SNI and NBS over FDR to be most effective at lower contrast-to-noise ratios (Zalesky et al., 2010).

The following takeaway points from Figure 6 pertain only to cases where the contrast forms a component and interconnects nodes that share a common property of interest (e.g., high centrality nodes): (i) as the CNR (μ_c) gets smaller, the better is the statistical resolution offered by NBS-SNI compared to NBS and FDR. (ii) The more pronounced is the node property of interest (parameter μ_nc) in the structural data, the greater is the statistical resolution offered by NBS-SNI. (iii) As the CNR is increased and μ_nc remains constant, the standard NBS and FDR will at some point offer a better resolution than NBS-SNI.

In sum, NBS-SNI offers a tighter control of the family-wise error rate than the other two methods in certain regimes by considering the nodes themselves, as opposed to considering only the edges. By giving importance to the nodes, the investigator can lower the parameter threshold t_edge to find significant abnormal connections of lower CNR. Conversely, the investigator could keep the same parameter threshold t_edge with NBS-SNI, thereby eliminating potential false positives that would result using NBS with the same value of t_edge. This choice would depend on the desirable trade-off between finding true positives and eliminating false ones. It is also important to note that if the structural connectivity matrices, or the anatomical property of interest are noisy, or if a contrast of interest does not span across representations, NBS-SNI might lower the statistical resolution offered by the standard NBS.

In Figure 7, we assessed the robustness of these conclusions by simulating structural connectivity matrices that might be more noisy, or of lower quality. We tested the case where NBS-SNI offered the greatest gain in statistical resolution in each row of Figure 6, that is, Figure 6C, D, and G.

For each pair of simulation parameters (CNR; K_nc(u, v)), we added three different levels of noise δ_extra-noise to the probability distribution P^structural(u, v) used to generate the structural connectivity matrices, that is, P^{noisy-structural}(u, v) = P^structural(u, v) + δ_extra-noise. The approach was the following: we computed P^structural(u, v) as before, but for each values P^structural(u, v) we added noise terms sampled from three different folded normal distributions (i.e., taking the absolute values of the normal distribution to prevent including negative entries in the probability distribution). The three folded normal distributions we tested are δ_extra-noise ∼ ∣𝒩(2.5, 2)∣, ∣𝒩(3.5, 2)∣, and ∣𝒩(4.5, 2)∣.

The simulation results in Figure 7 show that the gain in statistical resolution offered by NBS-SNI over NBS is most robust when the contrast-to-noise ratio (CNR) is low (Figure 7A), or alternatively, when the contrast on the central nodes (μ_nc) in the synthetic structural networks is high (Figure 7C). In the middle row, the CNR is too high and μ_nc is too low, and the standard NBS performs better than NBS-SNI for most of the discrimination thresholds t_edge. This example is a reminder for prospective user of the method that the extension we propose to the standard NBS is not intended as a replacement for the canonical method. Rather, we want to highlight the notion that in some regimes, NBS-SNI might offer a statistical gain over NBS. From these simulations, these regimes are characterized by a low CNR in the functional data, but also by an interplay of contrasts of interest in two different representations. In such cases, NBS-SNI can be used as a tool to leverage contrasts that span across representations.

The results presented in the ROC curves of Figures 6 and 7 offer a somewhat narrow view of the discriminating performances of the different methods, owing to the specific parameters employed in each simulation. A more comprehensive interpretation of these results is given in Figure 8, where the relative performances of NBS-SNI and NBS are shown when varying the number of samples and the overlap coefficient (OVL) between the Gaussian noise 𝒩(0, 1) and the contrast 𝒩(⁠$μe1$⁠, 1) (the standard deviation was again kept constant at σ = 1 for all realizations). The values shown on the x-axis represents the number of synthetic samples in each group. The OVL measures the similarity between two distributions via the overlapping area of their distribution functions. As for the ROC curves in Figure 6, N (N₁ = N₂) networks were generated in each group, with a contrast of 20 functional edges. Again, the position of the 264 nodes were taken from the from the UCLA Autism dataset (Rudie et al., 2013). Among these nodes, ∣n_c∣ = 5 were randomly selected to be the central nodes in the structural network.

To ensure a fair comparison of the true positive rate yielded by NBS-SNI and NBS, we controlled for the number of significant edges yielded by both methods. Their performances were compared based on subnetworks containing exactly 50 edges. To enforce this, the 200 (300) most significant t scores with NBS (NBS-SNI) were retained and the largest connected component identified. For NBS, the top 50 most significant t scores within the remaining connected component were retained. The same was done with NBS-SNI after taking the relaxed intersection (d = 1) with the central nodes (t_node = 0.5).

The orange circles represent a pair of parameters (number of samples, OVL) for which NBS-SNI offered a better performance (TPR_SNI > TPR_NBS), where their sizes and colors reflect the relative gain in statistical resolution compared to NBS. The relative importance (in both color and size) of each point is calculated as TPR_SNI − TPR_NBS. For visualization purposes, a colormap was not used for NBS, that is, when TPR_NBS > TPR_SNI (black triangles); only the size of the triangles reflects the relative gain in statistical resolution offered by NBS over NBS-SNI. For every point (circle or triangle) in Figure 8A, the simulation was repeated 500 times before the results were averaged to yield accurate visualization.

The results from Figure 8 show that with kind of synthetic data, the range of CNR for which NBS-SNI offers a statistical resolution gain over NBS decreases with the number of samples, and that the gain is most important when the CNR is low (i.e., when the overlap coefficient is high). Altogether, the results of Figures 6 and 8 indicate that NBS-SNI can improve the statistical resolution in certain regimes by considering two representations simultaneously. More precisely, we expect NBS-SNI to outperform NBS at identifying a contrast of relatively small contrast-to-noise ratio (CNR) in the functional network representations, where the nodes comprised in this contrast manifest themselves as high centrality nodes in the structural network representations of both groups. The intuition for proposing this kind of synthetic data is that some brain regions (nodes) could manifest as important information carrier (central nodes) in the structural representation of both the condition and the control group, while the arrangement and the strength of the connections emanating from these central nodes may differ in the functional representation. In other words, the synthetic data model employed is such that the identity of some high centrality nodes is common to both groups (condition and control) in the structural networks, but the functional connections between theses nodes may vary across the two groups. By relatively small, we mean that the CNR is smaller than the standard deviation of the modelled Gaussian noise.

Such a configuration may allow us to be less restrictive on the threshold applied to the set of edges in the functional networks. Since NBS-SNI is also concerned with central nodes, which are retained from the statistical test applied to the structural network representations, acknowledging central nodes in our framework might help us throw out some false positive edges that could have been retained from the functional edges’ two-sample statistical test. We posit that the false positive rate can be lowered in such a case with NBS-SNI, while the true positive rate can be increased as seen in Figure 8.

In this section, we present the results obtained with NBS-SNI applied to two real case-control studies. The first consists of individuals diagnosed with the autism spectrum disorder (ASD) and the second consists of individuals diagnosed with early psychosis (EP). All functional network representations are fully connected and weighted networks, where an edge represents the pairwise correlation between two brain regions. The procedure employed to yield predictions is detailed in the Methods. All the machine learning algorithms paired with NBS-SNI and NBS used to perform predictions are the ones implemented by the Python package scikit-learn (Pedregosa et al., 2011), without any hyperparameters optimization. In what follows, various machine learning (ML) algorithms were tested with the NBS-predict framework (Serin et al., 2021) for comparison purposes, but we want to emphasize that, depending on the research goal and the desired level of interpretability of the output networks, users should be mindful of which ML algorithms they base their findings on. Indeed, if interpretability of the outcome network is paramount to the research question, we recommend the investigators to utilize linear ML algorithms such as Linear Discriminant Analysis, Ridge Classifier, or Linear Support Vector Machine (SVM), as they allow for proper monitoring of the important edges (features) throughout the classification process.

This dataset consists of 29 individuals diagnosed with ASD and 19 healthy individuals (see Table 1 for the demographics). Each individual in both groups is represented by a structural connectivity matrix and a functional connectivity matrix (resting-state). Moreover, there exists a node correspondence between all 264 nodes (mapped according to the Power264 brain atlas of Power et al., 2011) across both groups and both representations. These structural networks are sparse weighted and undirected networks, where an edge represents the density of white-mater axonal fiber tracts between two brain regions. All the data was collected using a Siemens 3T Trio scanner. For the fMRI data, resting-state data was collected during 6 minutes, with TR = 3 s, TE = 28 ms and 3 × 3 × 4 mm voxels. The structural data was obtained from 8 minutes DTI scans with 2 × 2 × 2mm voxels. This publicly available dataset was downloaded from the USC Multimodal Connectivity Database (https://umcd.humanconnectomeproject.org).

We used NBS-SNI with d = 1 (i.e., edges at a maximal distance of 1 from the central nodes are retained), paired with the machine learning and cross-validation procedure described earlier to perform classification of brain networks in a case-control setup. Brain networks are classified in either the ASD or control group. The parameters K = 15 and r = 25 were employed for the repeated k-fold cross-validation.

The first node property investigated with NBS-SNI is the closeness centrality C_C (see Methods). The statistical test taking place on the measures of the nodes’ closeness centrality was the one-sample t test. The classification results obtained with different machine learning algorithms are presented in Table 2. The best prediction performance was yielded by the Logistic Regression algorithm (μ_performance = 0.670 ± 0.053) with the hyperconnected subnetwork, that is, a hyperconnected subnetwork in the autism group of brain networks, containing an average of k_hypo = 188 ± 20 functional connections across folds.

For the hyperconnected case with the Logistic regression algorithm, applying a weight threshold of 1 to the final adjacency matrix containing the average prediction weight of every edge across folds as its entries, 23 hyperconnected edges remained. A hyperconnected component of 11 edges was identified within these 23 edges. The hyperconnected component is depicted in Figure 9, alongside the list of its nodes and their associated degree. Interestingly, most of the connections are lateral (across hemispheres), with a larger proportion reaching the right putamen.

A variant of the closeness centrality, called the current flow closeness centrality (or sometimes information centrality) was also employed with NBS-SNI. As for the closeness centrality, the statistical test employed for this node property was the one-sample t test. The prediction performances with this node property are presented in Table 3. The best performance was yielded by the Ridge classifier algorithm (μ_performance = 0.734 ± 0.039) using the hypoconnected subnetwork, that is, a disconnected subnetwork in the autism group of brain networks containing an average of k_hyper = 227 ± 26 functional connections.

Applying a weight threshold of 0.99 to the final prediction adjacency matrix (obtained via the Ridge classifier algorithm) yielded 43 hypoconnected edges. Within these 43 edges, a hypoconnected components of 23 edges was identified, which is shown in Figure 10, alongside its nodes and their degree. Most connections seem to converge near the center of the sagittal plane of the brain. Notable is the large number of connections reaching the brain stem, the right cingulate gyrus posterior division, and the right precuneous cortex. Together, these three brain regions account for 46% of all the connections implicated between all 20 brain regions of the hypoconnected component.

Finally, the local information dimension was employed as the node property with NBS-SNI. Due to the nonnormal distribution of the computed nodes’ LID, a one-sample Kolmogorov–Smirnov test was used for this node property. To employ this node measure with the weighted networks of the structural representations, the sizes of the box l around a node was discretized into 15 equally spaced parts between the minimum distance and half of the maximum weighted distance between two nodes within the network. The prediction performances are presented in Table 4. The best performance was yielded by Decision tree classifier algorithm (μ_performance = 0.704 ± 0.063), using the hyperconnected subnetwork, containing on average k_hyper = 124 ± 18 functional connections.

Applying a weight threshold of 0.99 to the final prediction adjacency matrix yielded by the Decision tree classifier algorithm resulted in 17 hyperconnected edges. Within these 17 edges, the hyperconnected component of 10 edges was identified and is depicted in Figure 11, alongside the nodes and their degree. The majority of the nodes involved in the hyperconnected subnetwork are localized in the right hemisphere, connecting regions located solely in the frontal, temporal, and occipital parts of the brain.

Finally, the standard NBS was also employed to perform predictions. A range of parameter thresholds t_edge was tried, with t_edge = 2.5 being the one that yielded the best prediction performance. The classification results are presented in Table 5. The best prediction performance was yielded by the Gaussian naive Bayes algorithm μ_performance = 0.639 ± 0.026 with the hypoconnected subnetwork, containing on average k_hypo = 319 ± 33 edges.

To assess the significance of the improved prediction performance of NBS-SNI over NBS, permutations of individuals’ group memberships were performed. The prediction performances of both methods were calculated for 1,000 permutations. For each permutation, the group membership of every individual was randomly assigned before performing the nested cross-validation procedure. The p values were calculated as the number of times the maximum accuracy with random group memberships $μperformancenull$ was greater or equal to the empirical accuracy μ_performance, divided by the number of permutations (1,000). The results presented in Figure 12 show that the prediction performance of 74% yielded by NBS-SNI using the current flow closeness centrality was found to be statistically significant (p = 0.016) under the permutation test. For NBS, the accuracy of 64% was not deemed significant (p = 0.122).

The null distribution of prediction performances under group permutations is not centered at 50%, as one would expect. This indicates that prediction performances below ∼70% should be taken with a grain of salt with this dataset. This could be due to an artificial bias resulting from the small sample size of this dataset. Furthermore, this effect (i.e., a nonnegligible proportion of null prediction performances greater than 50%) is most pronounced with NBS in Figure 12B, therefore potentially showing a greater statistical discrimination offered by NBS-SNI with this dataset.

This dataset consists of 75 individuals diagnosed with nonaffective psychosis and 45 healthy control individuals (see Table 6 for the demographics) obtained from the Human Connectome Project (D. C. Van Essen et al., 2013). The fMRI data was parcellated using the Schaefer 200 7-networks parcellation (Schaefer et al., 2018). The anatomical measures of intrinsic curvature index and gray matter volume employed with NBS-SNI were calculated with FreeSurfer. The protocol scan sequences are T1w (MPRAGE) and T2w (SPACE) structural scans of 0.8-mm isotropic resolution. For the resting-state functional data, scans were obtained at 2-mm isotropic resolution, multiband (MB) acceleration factor of 8 and TR 720 ms, acquired twice: once with AP and once with PA phase encoding.

For this dataset, the anatomically derived measures of intrinsic curvature index and gray matter volume were employed as the node property within the NBS-SNI framework. In this case, the statistical test can be performed directly on the anatomical measure to probe for abnormal functional connections between certain anatomical regions. The parameters K = 15 and r = 10 were employed for the repeated k-fold cross-validation in the following results.

For the intrinsic curvature index, two independent one sample t tests were performed with the parameter p_node = 0.20, and the intersection between the two sets of nodes was retained. For the hypoconnected edges, the two sample t test parameter threshold was set to t_edge = 4.2. For the hyperconnected edges, this parameter had to be set to its minimal value t_edge = 0 in order to retain at least a few edges across all folds during the cross-validation process. This indicates that the early psychosis group had very few hyperconnected edges, while the control group had many more stronger functional connections. Once again, the parameter d = 1 was employed with NBS-SNI when taking the intersection between the nodes and the functional edges. The prediction performances are presented in Table 7. The best performance was yielded by Ridge classifier algorithm (μ_performance = 0.762 ± 0.015), using the hypoconnected subnetwork, containing on average k_hypo = 620 ± 200 functional connections.

After applying a weight threshold of 1, due to the large number of edges still remaining, only the top 10 nodes with the greatest degree were retained (corresponding to the 95th percentile) and the hypoconnected component of 16 edges depicted in Figure 13 was identified.

Using the gray matter volume as a node property, two independent one sample t tests were performed with the parameter p_node = 0.001, and the intersection between the two sets of nodes was retained. For the hypoconnected edges, the two sample t test parameter threshold was set to t_edge = 4.2. Once again, the parameter d = 1 was employed with NBS-SNI when taking the intersection between the nodes and the functional edges. The prediction performances are presented in Table 8. The best performance was yielded by Ridge classifier algorithm (μ_performance = 0.774 ± 0.029), using the hypoconnected subnetwork, containing on average k_hypo = 640 ± 191 functional connections.

After applying a weight threshold of 1, due to the large number of edges still remaining, only the top 11 nodes with the greatest degree were retained (corresponding to the 95th percentile) and the hypoconnected component of 20 edges depicted in Figure 14 was identified.

In the following results, a two-sample t test was applied to the node property of gray matter volume, rather than two independent one-sample t tests, as it was done for the previous results. Applying a two-sample t test on the node properties implies in this case that NBS-SNI probes for differences in gray matter volume across the condition and the control group, rather than searching for nodes that have high values of gray matter volume in both groups (i.e., resulting from taking the intersection of the suprathreshold nodes from both one-sample t tests). The reason for this choice was motivated by the large body of literature reporting loss of gray matter volume in individuals with schizophrenia (Gur et al., 1999, 2000; Hulshoff Pol et al., 2002; Vita et al., 2012). Applying a two-sample t tests on the node to investigate gray matter volume abnormalities between groups therefore seems appropriate.

NBS-SNI was employed with the parameters (⁠$tnodehyper$ = 0; $tedgehyper$ = 0) for the hyperconnected case and (⁠$tnodehypo$ = 1.2; $tedgehypo$ = 4.2). The prediction performances are presented in Table 9. The best performance was yielded by Ridge classifier algorithm (μ_performance = 0.776 ± 0.021), using the hypoconnected subnetwork, containing on average k_hypo = 903 ± 301 functional connections.

After applying a weight threshold of 1, due to the large number of edges still remaining, only the top 13 nodes with the greatest degree were retained (corresponding to the 95th percentile) and the hypoconnected component of 22 edges depicted in Figure 15 was identified. This hypoconnected subnetwork shows connections that were weaker (hypo), but also between nodes of lower gray matter volume in the early psychosis group. Using our framework, a thorough investigation of this condition could therefore investigate how functional connectivity is related with gray matter volume.

Finally, the standard NBS was also employed to perform predictions. A range of parameter thresholds t_edge was tried, with t_edge = 2.5 being the one that yielded the best prediction performance. The classification results are presented in Table 10. The best prediction performance was yielded by the Ridge classifier algorithm μ_performance = 0.727 ± 0.013 with the hypoconnected subnetwork, containing on average k_hypo = 9330 ± 990 edges.

After applying a weight threshold of 1, due to the large number of edges still remaining, only the top 10 nodes with the greatest degree were retained (corresponding to the 95th percentile) and the hypoconnected component of 42 edges depicted in Figure 16 was identified.

The results from the permutation test are presented in Figure 17. Overall, the prediction accuracies yielded by NBS-SNI with different anatomical properties, and the standard NBS were all deemed significant by this permutation test. When employed with gray matter volume as a node property, using a two-sample t test on the nodes’ statistical test, NBS-SNI offered a marginal gain of five percentage points of the standard NBS. In all cases, there was no overlap between the empirical prediction accuracy and its null distribution.

In this paper, a method called NBS-SNI was introduced, which is an extension to the Network-based statistic (NBS; Zalesky et al., 2010). By being supplemented with the structural network representation, or with anatomical measures, NBS-SNI can be used to identify abnormal functional connections and perform classification in case-control studies. More precisely, the framework that we presented sets to find the most important brain regions according to a certain node property of interest (e.g., centrality measures, anatomical measures), and probes for abnormal connectivity within or surrounding these thought-to-be important actors in the communication pathways of the brain. NBS-SNI is most well suited for cases when the contrast-to-noise ratio is relatively small in the contrast of interest. Furthermore, these connections forming the contrast should interconnect nodes that have some property that makes them stand out compared to other nodes; in that regime NBS-SNI is expected to offer a greater statistical resolution power than NBS by providing a tighter control of the family-wise error rate.

Using synthetic data, ROC curves were presented to show the regimes for which NBS-SNI is expected to offer a gain in statistical resolution over NBS at identifying a contrast. NBS-SNI was also applied on two real case-control datasets. The first one consisting of 29 individuals diagnosed with ASD and 19 healthy individuals containing both the structural and functional network representation of each individual. The second dataset comprised 75 individuals diagnosed with early psychosis and 45 healthy individuals consisting of resting-state functional networks, as well as the anatomically derived measures of intrinsic curvature index and gray matter volume.

In the fist dataset, using both network representations with NBS-SNI, a classification performance of 73% was obtained using the current flow closeness centrality of nodes, while the best performance obtained with NBS using only the functional representation was 64%. We assessed the validity of the prediction performances using permutation of the group memberships. For NBS-SNI, the prediction accuracy was deemed significant (p = 0.016), while the prediction accuracy of NBS was not deemed significant (p = 0.122).

For the second dataset, the anatomical measures of intrinsic curvature index and gray matter volume were used as a node property with our tool-box NBS-SNI. A classification performance of 77.6% was obtained with NBS-SNI using measures of gray matter volumes, while the best performance yielded by the NBS was 72.7%. All prediction performances were deemed statistically significant under the permutation test (p < 0.001).

We also want to emphasize the notion that NBS-SNI is a general framework, and its uses and formulation may not be limited to the ones presented in this paper. The choice of node properties to investigate is numerous. Using node properties such as centrality or anatomical measures, as it was done in this paper, seems most appropriate for large-scale neuroimaging data, since the volume of the parcellated brain regions does not vary much at this scale (Simpson et al., 2013), which is essential for the soundness of the subsequent statistical analysis undertaken with NBS-SNI, which relies on a close correspondence between nodes across individuals. Another avenue for the method could be to incorporate aggregated microscale information into this framework as the node properties. Interest and recent advances in multiscale brain imaging could provide researchers with a more comprehensive understanding of the structure-function relationships (Abdeladim et al., 2019; Desjardins et al., 2019; Schirner, McIntosh, Jirsa, Deco, & Ritter, 2018; Schulz et al., 2012). Indeed, such research endeavours should involve multiscale approaches, since the relationship between structure and function is believed to be heterogeneous across the brain, due to the different physiological and molecular microproperties of its different constituents (Khambhati, Sizemore, Betzel, & Bassett, 2018; Lariviere et al., 2019; Murphy et al., 2016; Suárez et al., 2020; van den Huevel & Yeo, 2017). Examples of microscale information that could be incorporated in a framework like NBS-SNI include gene expression, cytoarchitechture, neuron morphology, glial cells density/type, myelin content, etc. Animals such as rats and mice are well disposed to be subject to such multiscale neuroimaging studies (Desjardins et al., 2019; Schulz et al., 2012). These kinds of datasets will become more abundant in the coming years, and we believe NBS-SNI could be a useful tool to leverage them. Alternatively, if only one brain network representation is available, node properties could also be computed directly from the functional networks.

Besides, diagnosis of certain brain conditions, such as the autism spectrum disorder (ASD), can be challenging for clinical doctors, but also for the patients and their relatives. For example, a survey from 2016 in the UK shows that parents wait an average of 3.5 years between the time they approach clinical doctors to when they receive a diagnostic of autism for their children (Crane, Chester, Goddard, Henry, & Hill, 2016). Further, there are criticisms regarding the different diagnostic criteria for the different subtypes of the condition within the wider definition of the autism spectrum disorder (Tidmarsh & Volkmar, 2003). This is likely due to the overlap between the different subtypes of ASD; some are not well delineated, perhaps owing to the lack of a systematic quantification of the symptoms, brain morphology, and functional activity signatures implicated in the condition. We advance that network-based approaches applied in neuroscience, such as the method proposed in this work, are well positioned to address these issues, and could broaden the clinical doctors’ assessments with quantitative features.

Finally, another promising avenue for patients affected by a neurodevelopmental disorder who are refractory to conventional pharmacological treatment (Lindenmayer, 2000; Siskind, McCartney, Goldschlager, & Kisely, 2016) lies in alternative therapeutic approaches, such as transcranial magnetic stimulation (TMS) (Cash et al., 2021). This approach is becoming more precise in the targeting of the brain region(s) to stimulate. We posit that the identification of abnormal functional connections between influential actors in brain networks are promising candidates for customized targeted TMS treatments.

The formulation of NBS-SNI presented above only compares two properties between two groups. However, this could be extended to encapsulate more properties or more groups, or both. Multiple node or edge properties could be calculated for each of the two groups; the results of these could then be aggregated by Hotelling’s T² test into a test statistic that is distributed along an F-distribution (Al-Labadi, Fazeli Asl, & Lim, 2022). Alternatively, data from three or more groups could be compared using ANOVA (rather than the t tests used here), and this could then be used to select the nodes or edges for further comparison (note that this could also be used to compare more than two groups using traditional NBS). These avenues offer promise in larger scale future analyses, comparing more kinds of data from more varied datasets.

Furthermore, future analyses utilising NBS-SNI are likely to see its application in datasets with significantly larger sample sizes (e.g., the ENIGMA Consortium (Thompson et al., 2020), where studies contain thousands of participants). NBS-SNI is likely to offer statistical improvement over the canonical NBS in these datasets if neuroimaging data from multiple streams are used. However, one limitation of the present work is that the utility of NBS-SNI in large datasets has not been described. Due to the computational burden of generative network modelling (not the analysis of empirical data), we generated and evaluated the performance of NBS-SNI on 280 simulated networks. We therefore restricted our analyses of empirical datasets to those with a similar number of participants (given that some ground truth metrics had been established in this regime). We note that the performance of NBS-SNI worsened as the CNR and sample size increased (when compared to canonical NBS; Figures 6–8); further quantification of these effects would require large scale ground-truth simulation studies. As such, we believe that NBS-SNI is best suited to research paradigms where sample size is small (∼100–200 participants), CNR is relatively low, and there are complementary streams of neuroimaging data available. Although such paradigms remain common today, future usage of NBS-SNI in other regimes likely requires further modelling to validate its sensitivity and specificity.

It is important to acknowledge that it is difficult to provide definitive guidelines for the use of NBS-SNI at present. While we present an integrated framework for analyzing both nodes and edges, the choice of precisely how these properties should be analyzed remains with the investigator. These choices are present at different stages in the analysis pipeline: (i) when deciding on sources of information from which to generate node and edge properties (e.g., fMRI vs. diffusion MRI vs. structural MRI vs. alternative imaging modalities); (ii) when deciding which node or edge properties to use; (iii) when deciding what statistical thresholds to use to stratify significant nodes and edges; and (iv) when deciding which nodes and edges to retain (the parameter d).

Although no gold standard currently exists, some suggestions can be made regarding the options available in these domains. First, here we present examples using structural connectivity and anatomical measurements to define node properties. However, it would be reasonable to use metrics derived from alternative imaging modalities instead (e.g., as mentioned above). This provides practical benefits (given the cost and difficulty in acquiring tractographic results) and allows for the utilization of other streams of data which may be relevant to the condition under investigation. Second, the choice of which node or edge property to analyze may also be somewhat arbitrary. At this stage, we suggest that investigators consider the robustness of the imaging modality being used and the property being extracted from it. This should also be based on the extant literature regarding the condition under consideration, and multiple comparison corrections should be used if multiple node/edge properties are analyzed. Crucially, if one input matrix comes from a noisier source, NBS-SNI may end up losing statistical power. Third, as in much of the work that utilizes NBS and similar approaches, the choice of threshold may increase the likelihood of type I or II error (if the thresholds are too low or high, respectively). Since the original conceptualization of NBS, new techniques have been developed to allow for threshold-free analyses (e.g., Baggio et al., 2018). NBS-SNI could be augmented by these techniques in future work, but at present we suggest that investigators utilize thresholds derived from the effect size expected and the sample size available. Finally, we suggest utilizing the parameter d = 1 at present, which balances computational demands (given that computing shortest paths is not required) and exploration of subnetworks.

FN and AA acknowledge financial support from the Sentinelle Nord initiative of the Canada First Research Excellence Fund and from the Natural Sciences and Engineering Research Council of Canada (project 2019-05183).

Francis Normand: Conceptualization; Formal analysis; Investigation; Methodology; Writing – original draft; Writing – review & editing. Mehul Gajwani: Writing – review & editing. Daniel C. Côté: Methodology; Resources; Supervision. Antoine Allard: Conceptualization; Methodology; Project administration; Resources; Supervision; Validation; Writing – review & editing.

Antoine Allard, Sentinelle Nord, Université Laval (https://dx.doi.org/10.13039/100020862).