Global metrics.

Normalized global network metrics for each condition, across all thresholds, are reported in Figure 5. (As described in the Materials and Methods section, all the following network properties are computed on weighted matrices using weight-conserving algorithms. This approach, which is a departure from previous research in this field [25], [37], is motivated by the fact that binary matrices equally assigning a value of ‘1’ to all suprathreshold edges regardless of their connection strength are susceptible to false short paths which may significantly affect results [38] – see the Materials and Methods section for further discussion. In the Results and Discussion sections the ‘’ superscript is omitted for notation simplicity.)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Global metrics.

Average normalized global network metrics for each condition at each threshold (bars depict standard error). Abbreviations: Clustering, ; characteristic path length, ; small-worldness, ; efficiency, ; modularity, nQ. Significance level: ‘***’; ‘**’; ‘*’; ‘∼’; ‘ns’.

https://doi.org/10.1371/journal.pcbi.1003271.g005

Overall, the global repeated measures ANOVA indicated a significant effect of condition on normalized clustering (; , , ). The threshold factor was also significant (, , ), as expected, but did not interact with the condition factor (, ), indicating that the effect of the level of sedation on this measure is robust to the thresholding procedure. Follow-up pairwise comparison indicated that the effect of condition was mainly due to W and S exhibiting significantly less clustering than LOC and R, while no significant difference was found between either the first two or the latter two conditions. The follow-up 1-way ANOVAs (one per threshold) replicated the overall results at all thresholds, although it was only marginally significant at the lowest density threshold (i.e., ; see Figure 5, top row).

The level of sedation also affected the normalized characteristic path length (; , , ). Threshold level was also significant (, , ), but did not show any interaction with the level of sedation (). Pairwise comparisons indicated that the only significant difference occurred between LOC and all other conditions, with the state of unconsciousness exhibiting significantly greater (normalized) characteristic path length. These results were also replicated, at each threshold, with the 1-way ANOVAs (although the effect was only marginally significant at the lowest density threshold; i.e., ; see Figure 5, second row).

The effect of condition on small-world properties () generally mirrored that seen for , with a significant effect of condition (, , ), a significant effect of threshold (, , ) and no significant interaction between the two factors (, ). Follow-up pairwise comparisons indicated that, as for , the first two conditions (i.e., W and S) were not significantly different from each other, but both exhibited significantly smaller than the last two conditions (LOC and R), which were not significantly different from each other. In the 1-way ANOVAs, the same pattern was numerically detected at each threshold, but it was only significant at the three highest density thresholds (i.e., 21%, 26% and 31%), non-significant at the lowest density threshold (i.e., 11%) and marginally significant at the second lowest threshold (i.e., 16%; see Figure 5, third row).

Consistent with the results for , a significant effect of condition was also found on normalized global efficiency (; , , ). Threshold also exhibited the expected main effect (, , ), but again no interaction was observed with the level of sedation (, ). Pairwise comparisons indicated that in LOC efficiency is significantly decreased, as compared to all other conditions. The general pattern was replicated at each threshold individually, but it was only significant for the three highest density thresholds (i.e., 21%, 26% and 31%), marginally significant at the second lowest density (i.e., 16%; p = 0.06) and non-significant at the lowest density threshold (i.e., 11%; see Figure 5, fourth row).

Finally, normalized mean modularity (nQ) was also significantly affected by condition (, , ) and threshold (, , ), with no significant interaction between the two factors (, ). Pairwise comparison indicated that the effect of condition was mainly due to W and S exhibiting significantly less modularity than LOC and R, while no significant difference was found between either the first two or the latter two conditions. This same effect was seen in the follow-up ANOVAs at the three highest density thresholds (i.e., 21%, 26% and 31%), while only a marginal effect was seen at the second threshold (i.e., 16%) and no effect at the lowest threshold (i.e., 11%). With respect to the mean number of modules uncovered in each condition, at the global level we find a marginally significant effect of sedation (, ), although it exhibited a very small effect size (). Pairwise comparison indicated that the only difference was observed between S and LOC, with the latter showing a lower number of modules than the former. Consistent with the effect size statistic, however, although the same numerical trend was observed at each individual threshold (see in Table 2), it was never found to be significant in any of the follow-up 1-way ANOVAs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Mean number of modules for each condition at each density threshold (SD reported in parenthesis).

https://doi.org/10.1371/journal.pcbi.1003271.t002

Local metrics.

With respect to local metrics, results for nodal strength and local efficiency are depicted in Figure 6. The level of sedation had a significant effect on nodal strength across a wide variety of lateral and midline regions. Overall, two main patterns were detected. In some regions (shown in blue in Figure 6a), local strength was stronger (across thresholds and subjects) during W and R. This U-shape pattern was detected across a wide number of regions throughout the midline as well as in occipital, parietal and latero-ventral prefrontal-cortices (consistent with previous results [25]). In other regions (shown in yellow-red in Figure 6a), local strength was stronger during S and LOC, as compared to the other two conditions. This inverted U-shape pattern was mainly detected in temporal cortex, ventro-medial prefrontal cortex, and the apical aspect of prefrontal cortex. The effect of sedation on local efficiency () was, as expected (cf. [39]), consistent with the clustering results. All regions in which a significant effect of condition was detected exhibited greater efficiency (on average, across subjects and density thresholds) in LOC and R. As shown in Figure 6b, condition mainly affected the local efficiency of regions within the medial section of parietal and frontal cortices.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Local metrics.

Regions displaying a significant effect of condition on local metrics. (a) Nodal strength (yellow-red colors indicate regions in which degree was stronger, on average, for the S and LOC conditions, while blue-lightblue colors indicate regions in which degree was stronger for W and R). (b) Local efficiency (yellow-red colors indicate regions in which the measure is stronger, on average for the LOC and R conditions). Color intensity is assigned on the basis of the (FDR adjusted) p-value for the condition factor in the 2-way repeated measures ANOVA. (Surface rendering was performed using Caret [98].)

https://doi.org/10.1371/journal.pcbi.1003271.g006

Participants.

The whole dataset included twenty healthy right-handed volunteers (16 female) between the ages of 18 and 31 (M = 22.40, SD = 3.40). However, given the strong sensitivity of graph-theoretic measures to motion [56], [57], subjects presenting any displacement above 3 mm throughout any one scan (of the 4 each subject underwent) was dropped. In total, 8 subjects were excluded due to movement, reducing our sample size to 12 volunteers. For the remaining volunteers average absolute movement was well below half a millimiter (0.33, 0.34, 0.20, 0.26, for W, S, LOC and R, respectively) and did not significantly differ across condition (1-way repeated measures ANOVA; ; ). In addition, we remark that subject exclusion is, at present, a conservative strategy. Indeed, on the one hand, motion has been shown to affect graph theoretic measures [56], [57]. On the other hand, it has been shown that standard motion correction algorithms are insufficient corrections, and are liable to introduce spurious, but systematic, correlations [57]. Finally, novel more aggressive correction approaches censoring time-points [57], while correcting for motion and potential spurious correlations, are associated with the negative effects of reducing the number of time-points [65], [66] which has been shown to be associated with increases in the likelihood of high correlation, and to be a potential source of across condition or across subject bias [66]. As detailed below, in the following analysis we factor out motion as has been done previously [16], [25].

Conditions.

Volunteers underwent four resting-state fMRI scans. The across-scan variable was the level of consciousness, as clinically evaluated by the Ramsay scale [67]. The first scan was accomplished with participants being fully awake (wakefulness; W). In the second scan participants were sedated (S; Ramsay level 3) so that while their response to verbal command were slowed, they were still present. In the third scan participants experienced loss of consciousness (LOC; Ramsay levels 5–6), and exhibited no response to verbal instruction. Finally, the last scan was performed after participants had recovered (R; Ramsay level 2).

fMRI data acquisition.

Functional images were acquired on a 3 Tesla Siemens Allegra scanner (Siemens AG, Munich, Germany) with an Echo Planar Imaging sequence in 32 ascending slices (TR = 2,460 ms, TE = 40 ms, FOV = 220 mm, voxel size , and matrix size ). Because participants were acquired with a different number of volumes (varying between 196 and 350), correlation matrices were computed only on the first 196 volumes (8 min) for all conditions and subjects. For each participant one T1-weighted MP-RAGE image was also acquired (TR = 2,250 ms, TE = 2.99 ms, FOV = , , resolution isovoxel).

Preprocessing.

Functional and anatomical images were preprocessed according to the general procedures available in the 1000 Functional Connectome Project (http://fcon_1000.projects.nitrc.org/), and followed very closely the procedures employed in previous studies on the topic [16], [25]. First, the initial TRs of each functional dataset was removed. Second, data underwent slice-time correction, rigid-body adjustment for intra-run motion, brain extraction, 4 mm FWHM smoothing, band-pass filtering (), and removal of linear and quadratic trends. Nuisance signals, including motion parameters, white matter and CSF associated time-courses were partialled-out using a linear regression. Consistent with previous research [25], [37], global signal was not removed since it has been shown to be liable to introducing artifactual anti-correlations that can bias correlations differently in different parts of the brain depending on the underlying true interregional correlation structure, potentially introducing structure even where there is none [68], [69] and has been shown to suppress meaningful neural activity [70], [71], while regression of white matter and CSF signals has been shown to reduce many of the unwanted sources of noise that global signal regression is often used for [72]. Furthermore, with respect to graph theory analyses specifically, global signal regression has been reported to decrease reproducibility of both local and global topologically metrics [73]–[75]. The residuals of the regression were then co-registered to an MNI-space template (via 2-step registration using 6 degrees of freedom (dof) for within subject alignment of functional data to anatomical data and 12 dof to align the subject's anatomical image to the template). As part of this latter step, data were also resampled at isovoxel resolution.

Brain network constructions.

For each subject and each condition we constructed a graph representing a mathematical description of the brain as a functional network. A graph consists of a set of points and a set of lines connecting pairs of points [31]. In our framework, each point, also referred to as a node or vertex, corresponds to a specific brain region, or region of interest (ROI). Each line, also referred to as an edge or link, specifies the presence/absence of a connection between any two vertices and, for weighted graphs, the magnitude of the connection. (In our analysis the connection is the Pearson correlation statistic between each pair of nodes.) In what follows graphs are typically visualized in a familiar matrix heat-map presentation, where each row/column of the square matrix represents an ROI, and element (i,j) specifies the functional correlation between and .

Brain graphs were constructed in three steps. First, each individual data-set was parceled into 194 ROIs spanning cortex, sub-cortical nuclei, cerebellum and brainstem (see Figure 7). ROIs were defined independently, on the basis of a functional atlas that groups together spatially coherent voxels with homogeneous functional connectivity, at a desired resolution [76]. Specifically, we took the brain parcellation made available by Craddock and colleagues and employed that parcellation scheme to divide the brain in a set of 194 regions of interest. It is important to note that this procedure ensures that ROI selection is entirely independent of our data, thereby avoiding any form of bias in the analysis [77]. With the exception of atlas choice, this procedure follows exactly what has been done in previous research [25], [37], [78]. Choice of this parcellation scheme over previously used atlases such as the AAL [79] and the Harvard-Oxford Atlas [80], [81] is advantageous from several points of view. First, being functionally defined it clusters spatially proximal voxels by homogeneity of functional connections as opposed to clustering by anatomical position which, as exemplified by the case of the precentral gyrus ROIs in both the AAL and the Harvard-Oxford atlases, clusters together functionally very distinct subregions. Second, at our chosen level the Craddock ROIs have, collectively, almost twice the granularity as either structural atlas (i.e., 194 ROIs versus, 90 and 112 for the AAL and Harvard-Oxford atlases, respectively). Furthermore, when all atlases are equally resampled at 2 mm resolution, the ROIs obtained from the Craddock atlas appear much more homogeneous, in terms of size, than the other two (e.g., min to max voxel count range: 517, 5,005 and 8,076 for the Craddock, Harvard-Oxford and AAL atlases, respectively). After the functional parcellation of the brain into 194 ROIs, the average time-course associated with all voxels within each region was extracted. Finally, the average signal derived from each ROI was correlated to all other ROIs (using a Pearson's correlation coefficient; ) in order to obtain a square connectivity matrix (i.e., the graph). These connectivity matrices were the input for the following two analyses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. ROI selection.

Parcellation of brain data into 194 cortical, subcortical and cerebellar ROIs.

https://doi.org/10.1371/journal.pcbi.1003271.g007

Network classification.

Classifications were performed, on a pairwise basis, using the Pattern Recognition of Brain Image Data software package (PROBID; http://www.kcl.ac.uk/iop/depts/neuroimaging/research/imaginganalysis/Software/PROBID.aspx). PROBID employs a linear kernel Support Vector Machine (SVM) [82] algorithm to achieve classification; that is, it attempts to find a hyperplane that separates the matrices according to their class membership (e.g., W vs. S). This approach is particularly appealing in the domain of brain imaging because of its well-known benign properties in circumstances in which the number of input features (here, the number of elements of the adjacency matrix) exceeds the number of samples [83]. Use of a linear SVM algorithm has two important advantages. First, it reduces the probability of over-fitting the data (as compared to non-linear algorithms, since the linear kernel has only one parameter that controls the trade-off between having zero training errors and allowing mis-classifications; cf. [83], [84]); second, it allows direct extraction of the weight vector that defines the hyperplane, that is, a vector specifying which elements are most relevant for discriminating between the groups. In our context, the input to the classifier were the unthresholded correlation matrices, implying that ROI-to-ROI correlations were the features upon which classifications were performed. Following the requirements of the software, each subject's adjacency matrix was vectorized; conversely, the classification results were reshaped into the original matrix form. Classification performance was evaluated using an exhaustive “leave-one-out cross validation” scheme, as implemented in PROBID, resulting in mean accuracy, mean sensitivity, and mean specificity percentages for each classification. In this approach, N different training datasets are built by leaving out, in each, a different subject. The estimated SVM hyperparameter (i.e., the hyperplane best separating the two classes) from each training set is then separately validated against the one dataset that has been left out of the training. Finally, the results from each separate iteration of the training-validation procedure are averaged. Overall, it is important to stress that our approach avoids the most frequent pitfalls known to potentially bias machine learning approaches (see [83] for a detailed discussion). First, because feature extraction (i.e., computation of the adjacency matrices) was performed on a single subject basis, our subject-based cross validation strategy adheres to the independence assumption of training and validation set, thereby ensuring that (i) the estimate of the hyperparameter is indeed inferred solely on the properties of the training set, and has not been contaminated by the test set, and that (ii) the model evaluation is not biased by any information included in the training set. Second, because our rejection of artifacts (e.g., modeling of motion) is unconnected to class membership, there is no risk that the rejection schema might bias the classification results. Similarly, the decision to entirely reject a subject whenever one of the four sessions presented excessive motion also allows us to avoid the potential biases associated with unbalanced class frequencies. Finally, the use of an exhaustive leave-one-out validation strategy provides an almost unbiased estimate of the generalization error [83], albeit at the expense of an increased variance of the estimator, thereby avoiding the possibility of excessively optimistic estimates of the generalization error. Overall, five classifications were computed. First, we attempted to classify the two states of wakefulness (W & R) versus the two states of altered consciousness (S & LOC). Then, we performed three pairwise classifications aimed at distinguishing contiguous states; namely, W vs. S; S vs. LOC; LOC vs. R. Finally, we compared the two states of wakefulness to each other (i.e., W vs. R). The statistical significance of each classification accuracy was assessed with a non-parametric approach based on repeated label permutation (10,000 times). The p-value from the permutation was derived by counting the number of times the true (that is, non-permuted) accuracy was greater than the accuracy derived after permuting the labels and dividing by 10,000. Familywise error rate was controlled for with a Bonferroni correction.

To visualize the results of the SVM classification for each comparison, we created “importance maps” (also referred to as “evidence maps”) [85], [86]. These maps, depicted in Figure 3, were obtained in two steps. First, we separated the elements of the weighting matrix w, which defines the hyperplane best separating two classes (e.g., LOC versus S), into two sparse matrices. The first sparse matrix, , included all the elements contributing to classifying the first condition (e.g., LOC); the second sparse matrix,, included all elements contributing to classifying the second condition (e.g., S). (It should be noted that while there is no absolute scale for quantitatively interpreting the importance maps, the greater the value for a given link, the greater its influence on the classification.) Second, we multiplied each sparse matrix by the average adjacency matrix across all subjects correctly classified. The resulting sparse matrices, which are depicted in Figure 3, thus show the average value, across all subjects correctly classified, of each edge contributing to the classification of a given condition scaled by how much it contributed. Positive values therefore indicate positive correlations (i.e., edges) contributing to a classification, scaled by how much they contributed to the classification; similarly, negative values indicate negative correlations (i.e., edges) contributing to a classification, scaled by how much they contributed to the classification.

In parallel to the SVM algorithm, all classifications were also computed using PROBID's Gaussian Process Classifier (GPC) [87] algorithm, which yielded nearly identical results.

Finally, to assess whether the pattern of connections which allowed classifying the transition from wakefulness to sedation (i.e., S vs. W) differed from the pattern of connections which allowed correctly classifying the transition from sedation to loss of consciousness (i.e., LOC vs. S), we performed a test comparing the distribution of classifying connections in the two comparisons. More in detail, for each of the two classifications, we divided the top 1% of edges contributing to correct classification into 15 kinds of connections (e.g., cortico-cortical (positive), cortico-cortical (negative), thalamo-cortical, striato-cortical, etc; see Tables S1 and S2). Counts for the two conditions were then entered into a contingency table and thus submitted to a test comparing the two distributions.

Network analysis.

To characterize local and global properties of the brain networks during normal wake, sedation, loss of consciousness and recovery, we employed a set of metrics derived from graph-theory [27], [54], [78]. Prior to the analysis, connectivity matrices were thresholded to avoid the possibility that weak and non-significant links representing spurious connections may obscure the topology of strong and significant connections [28]. Thresholding was performed on a proportional basis thus retaining, for each individual matrix, only the strongest t% (positive) edges. To date there is no agreed algorithm to select a unique threshold at which to perform graph-analysis, hence we followed one of the many procedures employed previously. Specifically, the lower boundary was selected to ensure the averaged degree was not smaller than , where N is the total number of nodes in the graph (i.e., N = 194). This lower boundary guaranteed that the resulting networks were estimable networks [30]. The upper boundary was selected to ensure that mean small-worldness (see below for the definition) in each condition was not smaller than 1.0. These two constraints yield an interval of density thresholds between 11% and 31%, which we sampled in steps of 5% (resulting in the following five density thresholds: 11%, 16%, 21%, 26%, 31%). The issue of thresholding is addressed in more detail in the Discussion.

Several global and local metrics were computed for each connectivity matrix, employing the Brain Connectivity Toolbox (BCT) [28]. First, we assessed large-scale network functional integration of brain networks, which captures the extent to which information from distributed regions can be rapidly combined, by measuring their characteristic path length and global efficiency. The characteristic path length (L) of a network is defined as the average length of the path uniting each pair of nodes within the network [30]. Since our networks are based on functional data, paths can be conceived as sequences of statistical associations (i.e., correlations) between nodes rather than physical existing/non-existing paths as might be the case with structural measures. The most common and mainstream interpretation of characteristic path length in networks is as a measure for global integration of information between topologically distant brain regions (see Discussion for more details). To date, the vast majority of graph analysis of functional data have been performed on binary graphs, that is, connectivity matrices in which all elements that survive thresholding are set to 1, regardless of their original value, and all other elements are set to 0. In this procedure, however, once matrices are binarized weak edges exert the same influence on network measures than strong edges. In the context of path length, for example, where strongly correlated nodes can be intuitively interpreted as (functionally) “closer,” binarization might result in erroneous creation of long-distance shortcuts, which can obscure the architecture of strong connections and, ultimately, important across-group differences [38]. Thus, following the development of metrics that make use of weight information [88], [89], we employ “weighted” measures (denoted, in the following formulae by the superscript ‘’). Characteristic path length was thus computed as follows:where represents the set of all nodes within a network, the number of nodes, and is the shortest path uniting nodes . The average inverse of L is thus a measure of global efficiency () of a network [39]. A different aspect of a network is the degree to which it exhibits functional segregation, which captures the extent to which subgroups of nodes densely interconnected with each other can carry out specialized processing. To quantify the presence of clusters of regions (e.g., modules) with strong statistical dependencies suggestive of segregated neural processing we computed the clustering coefficient (C), a measure of the ‘cliquishness’ of networks' subdivisions defined as the average fraction of a node's neighbors that are also neighbors of each other [30]. In its weighted version, clustering is defined as [89]:where is the degree of a node (i.e., the sum of the weights of all nodes connecting to node ), and represents the geometric mean of all triangles (i.e., triads of reciprocally connected nodes) surrounding node . Similarly, modularity (Q) measures the extent to which a network can be optimally subdivided in a number of non-overlapping groups of regions which have high within-group connectivity and low across-group connectivity [90]. As for above, we employed the weighted version of the algorithm [88]:where is the weight of the connection (i.e., correlation) between nodes and , designates the module within which node is contained, and is equal to if nodes and are part of the same module , and otherwise. To address the issue of degeneracy, we followed a three step procedure previously suggested [28]. Specifically, computation of modularity was refined by first applying a finetuning algorithm to the above calculation, and then running a final probabilistic tuning algorithm on the finetuned result. In addition, because the calculation of modularity is heuristic and variable across iterations, we performed the three steps above (modularity, fine-tuning and probabilistic-tuning) 50 times for each matrix, and then employed the mean modularity value across the 50 iterations for subsequent analyses (we note that we also performed the analyses on the median modularity, across the 50 iterations, and it yield qualitatively identical results).

It is important to note that it is often not possible to interpret the significance of the absolute values of the above metrics. Thus, the metrics derived for each network are typically normalized by the average metric derived from 100 random networks (with matched size and degree distribution). The normalized characteristic path length () is thus the ratio of the characteristic path length () of a brain network and the average characteristic path length of 100 random networks (). Normalized clustering coefficient (), global efficiency () and modularity (nQ) are defined analogously. In addition, the ratio of normalized clustering and normalized characteristic path length is often referred to as the small world property () of a network. Small-world networks are networks that are neither completely random nor completely regular; rather, they lie in-between these two extremes and exhibit a high degree of clustering, like regular networks, and a low characteristic path length, like random networks [30]. Several studies to date have observed small-world characteristics in the human brain [26], [91].

In addition to the metrics above, we also characterized the brain networks on a node-wise basis. First, we computed the strength () of each node, a basic measure of node centrality defined as:Nodes with a high strength are nodes that are (functionally) close to many other nodes in the graph. In general, measures of centrality are based on the intuition that nodes participating in several short paths within a network are important controls of information flow [31]. Node strength can therefore be interpreted as an index of communication activity within a network. Finally, we also calculated the local efficiency () of each node to get a local measure of segregation. In its weighted version, local efficiency for node i is defined as [28], [39]:where is the length of the shortest (weighted) path between and that contains only neighbors of . This metric can be interpreted as measuring how “fault tolerant” the system is since it captures how efficient the communication is between the first neighbors of when is removed [39]. The present extension of the local efficiency to weighted networks differentiates the influence on presumed fault tolerance of strongly connected and weakly connected neighbors, and broadly parallels a generalization of the clustering coefficient to weighted networks [28], [89].

Statistical analyses were performed in Stata (SE, Version 13, Stata Corp) and Matlab (Mathworks Matlab, version R2013a). The effect of condition on each metric was assessed using a repeated measures ANOVA model with condition (4 levels: W, S, LOC and R) and threshold (5 levels: 11%, 16%, 21%, 26%, 31%) as factors. Following convention [92], where the Machuly's test indicated violation of the sphericity assumption, degrees of freedom were calculated with the Huynh-Feldt correction for moderate violations (; denoted with the subscript) and the Greenhouse-Geisser for more severe violations (; denoted with the subscript). Because of the relatively small size of our sample, we also report the partial omega-squared (), a measure of effect size in the population that is resistant to sample size and therefore provides a better measure of the magnitude of the effect between conditions [93], [94]. Unlike other metrics of effect size that are significantly upwards biased (e.g., ), does not results in inflated estimates even for small sample sizes [93], [95]. This statistic is usually interpreted as the percent of the dependent variable's variance accounted for by the effect in the population with other non-error sources of variance being partialled out, and is considered to indicate, for values of 0.01, 0.06, 0.14, small, medium and large effects, respectively [96]. Pairwise comparisons among the different conditions were performed on the basis of estimated means using a Sidak adjustment for multiple comparisons. The 2-way repeated measures ANOVA was followed-up with a set of 1-way repeated measures ANOVA, one per each threshold, with condition as the only factor. For global measures, significance was assessed against a conventional 0.05 p-value criterion. For local measures, due to the multiple comparisons issue, we employed an FDR-adjusted criterion of (following [97]).