Temporal segmentation of EEG based on functional connectivity network structure

In the study of brain functional connectivity networks, it is assumed that a network is built from a data window in which activity is stationary. However, brain activity is non-stationary over sufficiently large time periods. Addressing the analysis electroencephalograph (EEG) data, we propose a data segmentation method based on functional connectivity network structure. The goal of segmentation is to ensure that within a window of analysis, there is similar network structure. We designed an intuitive and flexible graph distance measure to quantify the difference in network structure between two analysis windows. This measure is modular: a variety of node importance indices can be plugged into it. We use a reference window versus sliding window comparison approach to detect changes, as indicated by outliers in the distribution of graph distance values. Performance of our segmentation method was tested in simulated EEG data and real EEG data from a drone piloting experiment (using correlation or phase-locking value as the functional connectivity strength metric). We compared our method under various node importance measures and against matrix-based dissimilarity metrics that use singular value decomposition on the connectivity matrix. The results show the graph distance approach worked better than matrix-based approaches; graph distance based on partial node centrality was most sensitive to network structural changes, especially when connectivity matrix values change little. The proposed method provides EEG data segmentation tailored for detecting changes in terms of functional connectivity networks. Our study provides a new perspective on EEG segmentation, one that is based on functional connectivity network structure differences.


Proposed algorithm
Functional connectivity network construction Given a segment of EEG data, the network connectivity matrix first needs to be calculated (Fig. 2A).Many connectivity metrics have been proposed, and among them commonly-used ones include: Pearson correlation coefficient (PC) 24 , mutual information 25,26 , phase-locking value (PLV) 27,28 , and phase lag index (PLI) 28 .We here use PC and PLV, though our segmentation method is agnostic to the particular choice of connectivity metric, and the merits of choices are not the focus of our study.A functional connectivity network is then constructed based on the connectivity matrix.The nodes represent brain regions surrounding corresponding EEG electrodes and edges represent the strength of communication between brain regions.When using a threshold on connection strength, the connectivity matrix becomes an adjacency matrix and the network becomes an unweighted network.A weighted network can be constructed if not using a threshold and keeping the original connection strengths.We work with undirected weighted networks in this study, so as to retain more information.

Computing distance between graphs
Graph distance measures (or graph similarity) describe the dissimilarity (similarity) between networks, quantifying topological differences between networks.Graph distance measures have long been a focus in network science and are widely used in many fields, such as pattern recognition 29 , model selection 30 , network classification and clustering 31 , anomaly, and discontinuity detection 32 .Previously proposed measures include [33][34][35][36][37] .We here build an intuitive and modular framework for graph distance by calculating a vector to summarize each graph and then computing distances between vectors representing different graphs.The vectors contain a value per node, where the value is an index of the node's importance.This framework allows different node importance indices to be plugged in, resulting in different graph distance measures (abbreviations with the DM prefix).
Given two time periods T 1 and T 2 and EEG data with n channels, connectivity matrices W(T 1 ) and W(T 2 ) are calculated from the time periods, respectively.Weighted networks G(T 1 ) and G(T 2 ) are constructed from the connectivity matrices.The indices of node importance in networks G(T 1 ) and G(T 2 ) are represented by n-length vectors c(T 1 ) = (c 1 (T 1 ), c 2 (T 1 ), . . ., c n (T 1 )) and c(T 2 ) = (c 1 (T 2 ), c 2 (T 2 ), . . ., c n (T 2 )) .We calculate graph distance d(c(T 1 ), c(T 2 )) using the Euclidean distance between two node importance vectors, as follows: We here consider six indices of node importance for weighted networks 38 , summarized in Table 1.DC, CC, EC, BC, CCO, and NNC denote degree centrality, closeness centrality, eigenvector centrality, betweenness centrality, clustering coefficient, and nearest neighbor centrality, respectively.Applying Eq. (1) to each node importance index gives each of the six corresponding graph distance measures: DMDC, DMCC, DMEC, DMBC, DMCCO, and DMNNC, respectively.These six indices are briefly compared in the "Discussion".A schematic of the procedure for calculating graph distance is shown in Figure 2C.

Sliding window change detection
To perform temporal segmentation of multi-channel EEG data, we designed a sliding window change detection method based on a previously published method 23 that compares a sliding window (representing new data) with a growing reference window (Fig. 3, for differences from their algorithm, see the "Discussion").The parameters of the algorithm are listed in Table 2 (for how these are set, see "Method parameters and evaluation", a brief discussion of the two most important parameters is in the "Discussion").In brief, the algorithm calculates the distance between the network in the reference window and that in the sliding window; if this distance is abnormally large (an outlier among samples of distance values), a boundary is declared.The algorithm steps are as follows: (1) Table 1.Definitions of node importance indices.Abbreviations in parentheses are names given to distance measures based on these node importance indices.W is the connectivity matrix of the network with elements w ij .n is the number of nodes.N i is the set of all nodes connected to node i. N is the set of all nodes in the network.k i is the number of the set N i .d ij is the length of the shortest path between node i and j. g kj is the number of shortest paths from node k to node j, and g kij is the number of shortest paths that pass through node i from node k to node j. x i is the eigenvector of the max eigenvalue of matrix W.

Index Definition
Degree centrality (DMDC) Closeness centrality (DMCC) Eigenvector centrality (DMEC) Betweenness centrality (DMBC) Clustering coefficient (DMCCO) The time periods of the reference window and the sliding window are denoted as T r and T s , respectively.(ii) EEG data in the reference window is placed in X r .EEG data in the sliding window is placed in X s .(iii) Functional connectivity matrices of X r and X s are calculated (with the chosen connectivity measure, e.g.correlation) and stored in A r and A s , respectively.(iv) Undirected weighted networks G r and G s are constructed from A r and A s , respectively.(v) Node importance values for each node in G r and G s are calculated with the chosen node importance index (e.g.closeness centrality), giving c(T r ) and c(T s ) , respectively.(vi) Graph distance measure d(k) is calculated by equation (1) with d(c(T r ), c(T s )) , i.e. the L2 distance between the node importance value vectors.(vii) Since we detect change by detecting outliers, a sufficient sample size is necessary to obtain reliable distributions for outlier detection.Hence, we set this criterion: If k is less than W KDE d , then (not enough data): k is increased by 1 and d(k) is appended to list L; t i + W − W v is the time corresponding to this d(k) entry in L; the reference window grows by W p , i.e. its length becomes W = W + W p , while its starting time does not change; the sliding window slides forward in time by W p ; go to step ii.(viii) Outliers of list L are calculated (see "Outlier detection").(ix) To declare a boundary, an outlier must have been found.k must also be at least the required number of comparisons W d .Lastly, the distance d(k) to the current sliding window should be no greater than the detected outlier in L. Hence, we check for these three criteria: (1) outliers in If they are all true, we detect a boundary (the largest outlier's corresponding time is the detected boundary time) and go to step i. (x) d(k) is appended to list L; the t i + W − W v is the time corresponding to this d(k) entry in L; k is increased by 1; the reference window grows by W p ; the sliding window slides forward in time by W p ; go to step ii.

Outlier detection
Outlier detection is a key step in the segmentation algorithm: it finds time points with larger than expected graph distances between windows.Considering that the values in list L may not be from a normal distribution, we use a non-parametric method for testing outliers.First, we use kernel density estimation (KDE) to obtain a density function.Then a cumulative distribution function P(d) is estimated.Let P(d < d threshold ) = P KDE , and the threshold d threshold can then be obtained.P KDE is a parameter that specifies the percentile threshold for outliers.All elements of list L greater than d threshold are outliers.
To verify the necessity of using a non-parametric method to detect outliers, we performed normal distribution tests for L (computed using GD) on real EEG data from 28 participants.Every time outlier detection was required, a normal distribution test was performed.The null hypothesis of normal distribution (significance level 0.05) was rejected on average (across participants) 35.3% of the time.

Comparison algorithms
To compare against our proposed framework, we plug in (dis)similarity measures proposed in previous work 14,23 .While their original algorithms work on raw data matrices, we modify them to work on functional connectivity matrices in the interest of fairness, as our simulations focus on functional connectivity changes.Below we present the connectivity-matrix-based (dis)similarity measures used for comparison.In contrast to our proposed method, these approaches focus more on the difference in connection values, rather than graph structure.The EEG data matrix based methods for segmentation 14,23 use the Kolmogorov-Smirnov (KS) test and the Grassmann distance.The KS test and Grassmann distance are used here on two functional connectivity matrices, instead of two data matrices.
For both the KS and Grassmann approaches, the connectivity matrices W(T 1 ) and W(T 2 ) and their eigen- decomposition are calculated to determine their column space and weights of corresponding eigenvectors.The significant eigenvectors are selected by the elbow method to construct feature subspaces F(T 1 ) and F(T 2 ) , respectively.The Grassmann distance between subspaces is calculated as follows: where k = min(dim(F(T 1 )), dim(F(T 2 ))) , dim(•) is the dimension of a subspace.θ i (F(T 1 ), F(T 2 )) denotes the ith principal angle between subspaces.For the KS test 23 , each column of W(T 1 ) and W(T 2 ) is projected onto the feature subspace F(T 1 ) of the reference window, and the residual of the projections are gathered into two error matrices denoted e(T 1 ) and e(T 2 ) , respectively.Then, the Kolmogorov-Smirnov test is used to test the difference between e(T 1 ) and e(T 2 ) , and the final similarity measure is the test's p-value.A flowchart of the matrix comparison process is shown in Fig. 2B.
In the sliding-window change-detection algorithm, the changes needed for KS and Grassmann approaches are as follows.Grassmann distance: In steps iv-vi, feature subspaces of the matrices A r and A s are obtained and Grassmann distance d(k) is calculated by Eq. ( 2).Kolmogorov-Smirnov (KS) test: In steps iv-vi, matrices A r and A s are projected onto the feature subspace of A r .The residuals after projection are denoted e(T r ) and e(T s ) , respectively.The p-value of the KS test between e r and e s is stored in d(k).In outlier detection for the KS test, the smaller the p-value, the greater the difference between two matrices.Thus, in step ix: the condition is reversed,

Simulation results
We first test our segmentation method on simulated EEG data where ground truth of the dynamics of functional connectivity is known.The purpose of this test is to verify that the segmentation method finds the time points where connectivity structure changed.We also vary the difficulty (i.e.magnitude of changes) and examine the effects on accuracy of finding change points.We compare our proposed method with that based on connectivity matrix (edge weight) differences.

Simulation design
In this section, our segmentation method is evaluated on simulated EEG data.We simulated EEG data with N = 32 channels, T = 60 s time length, fs = 100 Hz sampling frequency, and a total of S = 200 repetitions.The EEG simulation process consisted of first generating several different functional connectivity networks (so that switching among them gives dynamic network structure) , i.e. specifying target connectivity matrices and then generating time-varying EEG data with the target connectivity matrices.We simulated two kinds of changes in network structure, as follows: Simulation 1: dynamic network with changing communities In this simulation, the number of communities changed from 2 to 3 and back.In the interval [1, 20] seconds, the functional connectivity network had two community structures.Community 1 consisted of nodes 1-10 and Community 2 consisted of nodes 15-32.The strength of the connections within the communities was higher (drawn from a normal distribution N(0.2, 0.1)) than that between the communities (drawn from N(0, 0.05)).To add random diversity, 22 edges and 76 edges were selected from Community 1 and Community 2, respectively, and random weights (N(0.2, 0.1)) were added to the original edge weights (the R-MAT algorithm 39 , one of the most commonly-used network generation models; it models graph structure by a degree distribution that follows a power-law distribution, approximating the properties of real-world networks; R-MAT has been previously applied in simulating dynamic brain networks 14 .).The network in the interval [40, 60] s is the same as that in the interval [1, 20] s.
In the interval [20, 40] s, the functional connectivity network had three community structures.Community 1 and Community 2 were the same as in the interval [1, 20] s.Community 3 consisted of nodes 11-14.Weights of the edges within Community 3 were drawn from N(0.2, 0.1), and the edges connecting to other communities were drawn from N(0, 0.05).According to the R-MAT algorithm, 4 edges in Community 3 were selected to add random weights ( N(0.2k s , 0.1) ), where the signal strength parameter k s ranged from [0, 1].There is a proportional relationship between the k s value and the signal-to-noise ratio (SNR), as shown in Table 3.The calculation of the SNR is based on the connectivity matrix using a previously published method 14 .The larger k s , the greater the change of the connectivity matrix, and the greater the SNR.The configuration of the connectivity matrix and its corresponding network are illustrated in Fig. 4A.

Simulation 2: dynamic network with central hub
Central hubs are a key feature of human brain structural organization, play an important role in function 40 , are involved in changes in brain networks 41,42 , may be useful in decoding brain activity 43 , and may dynamically alternate 44 .Considering their importance, we simulated EEG data where a central hub connecting communities appeared and disappeared.
In the interval [1, 20] s, the connectivity network was the same as that in the interval [1, 20] s of the previous simulation.The difference from Simulation 1 is that Community 2 is set to nodes 13-32.The number of additional weighted edges within Community 2 was changed from 76 to 95. Between [20, 40] s, the network had three communities: Community 1 consisted of nodes 1-10, Community 2 consisted of nodes 13-32, and Community 3 consisted of nodes 10-13.Note that the overlapping nodes make Community 3 a central hub connecting Community 1 and Community 2. The weights of Community 1 and Community 2 were the same as in the interval [1, 20] s.The weights within Community 3 were drawn from N(0.2, 0.1), and the edges connected to other communities were drawn from N(0, 0.05).Four edges within Community 3 were selected and random weights were added using the R-MAT algorithm.The randomly added weights were drawn from N(0.2k s , 0.1) , where the range of the signal strength parameter k s was [0, 1].The configuration of the connectivity matrix and its corresponding network are illustrated in Fig. 4B.

Time-varying EEG data generation
There are many methods to simulate time-varying EEG data.For example, the SEED-G toolbox 45 can simulate EEG data with preset phase-locking values (PLV), phase lag index (PLI), and directional phase lag index (dPLI).However, EEG data with correlation between more than 3 channels cannot be generated due to high 2, where a central hub appears (as shown) and subsequently disappears.Note that the central hub differs from a regular community in that it has one node which is well-connected to each community.
computational memory load.The SimMEEG toolbox 46 can set the number of network nodes (channels), connection strength, network density, and other parameters, but it cannot simulate EEG data with community structure.Šverko 47 proposed a simulation method that can generate EEG data with community structure, but the community structure is not stable enough over time and connections between communities cannot be preset.Moiseev 48 proposed a simulation method that allows specified Pearson correlation between channels.Considering computational load and stability of the resulting correlation matrices, we use this method.In this method of generating EEG signals, time courses are generated to have user-specified mutual correlations, and generated signals consist of a combination of evoked signals and (randomized) oscillatory signals.To meet the method's requirement for positive definite connectivity matrices, we use an iterative spectral method to bend non-positive-definite symmetric matrices to be positive-definite 49 .The bending has little effect on the connectivity matrix weights (Supplementary Figure S1).
To make the properties of the Pearson correlation connectivity matrix closer to those of real EEG data, the connectivity matrices generated in the previous sections need to be adjusted.These adjustments do not change the network structure.First, we normalize: given a connectivity matrix W of size n × n , the connectivity matrix is normalized so all values are [0, 1]: Second, we adjust the magnitude of weights, since in real data, Pearson correlation coefficients are not particularly large.We reduce the correlation coefficients in the simulations by multiplying by 0.6 50 , resulting in values in [0, 0.6].Third, the connectivity matrix is symmetrized: Fourth, the diagonals of the matrix are changed to 1:

Method parameters and evaluation
The parameters of the segmentation algorithm need to be set to the same values across compared algorithms to ensure fair comparison.We first determined the sensitivity of parameters (i.e. which parameters have a large impact on segmentation results) and found that W r and P KDE have a relatively large influence on segmentation performance.Thus, it was necessary to perform a grid search for the optimal values of these two parameters (Supplementary Figure S2).For the relatively insensitive parameters, we set W s = 2 s, W v = 1 s, and To evaluate the performance of each distance measure when combined with the segmentation method, we use several metrics 23 , as follows.Success rate ( p succ ): ratio of number of successfully detected boundaries to the total number of ground truth boundaries.As the length of the sliding window was 2 s, if the detection boundary was within 1 s of the real boundary distance, it was considered a successful detection.Failure rate ( p fail ): ratio of number of falsely detected boundaries to the total number of ground truth boundaries (s).Aggregate rate ( p aggr ): the difference between success rate and failure rate.Average estimated displacement ( µ ED ): mean of the absolute distance between the detected boundary and nearest ground truth boundary.Here we include all declared detections.The standard deviation of estimation displacement ( σ ED ): variability in the displacement around the average displacement µ ED .

Segmentation results
Eight measures (six choices of node importance indices describing graph structure for our proposed method, two choices of matrix difference measures representing previous work) were plugged into the segmentation framework and tested on the simulated EEG data.Among these, GD and KS are matrix-based measures which represent the traditional approach to segmentation.We also evaluated the Source-Informed Segmentation (SIS) method 23 in preliminary analysis, but performance on our simulated data was poor (results not shown).
The aggregate rate, the success rate, the failure rate, and the average and standard deviation of the boundary displacement at different k s values are shown in Figs. 5 and 6.From the overall trends we can see that, with an increase of signal strength k s (the "magnitude" of changes), the performance of distance measures other than KS gradually improved.In Simulation 1, DMCC and DMDC substantially outperformed GD and KS, the two traditional, matrix-based measures.This indicates that with these two node importance indices, networkstructure based segmentation (using indices related to graph structure) performs better than functional connectivity matrix based segmentation (directly comparing the edge weights).
When k s was low, KS and DMCCO also performed well, but segmentation performance did not increase much with k s .In Simulation 2, DMCC and DMBC were better than GD and KS, especially when signal strength k s was low.DMDC performed slightly better than GD.This again indicates network-structure based segmentation performs better than matrix based segmentation.KS performed poorly in this simulation.The curves for DMEC and GD were found to overlap for both simulations.This is because DMEC and GD are both calculated based on eigenvectors.
To give specific numbers to make the results more intuitive, in Simulation 1, at SNR of about 3dB (signal about twice as large as noise), the DMCC method had a success rate of finding change points of around 88%, with mean www.nature.com/scientificreports/Similarly, let s t i ,t i+1 be the number of p t i ,j larger than p t i .We define a difference ratio as follows: The higher the p diff value, the better the segmentation performance.The calculation process of p diff is shown in Fig. 8A.
Since the easy condition and difficult condition were performed in two blocks with a rest in between, we segmented data from the two conditions separately.In the experiment, there were sound stimuli (SS) at random times, which produced event-related potentials in the EEG data.Therefore, we also used the sound stimulation times as segment boundaries for comparison.We also added a naive segmentation method for comparison, where each segment was the same length (SL), 800 samples (for a similar total number of segments as SS).
We repeated segmentation with Pearson correlation (PC) and phase locking value (PLV) as two choices for the connectivity matrix calculation method.The parameters of the segmentation algorithm were set as: The parameters (W r , p KDE ) had a relatively large effect on p diff and the number of segments found.For these two parameters, we did a grid search, and the results for GD and DMDC are shown in Supplementary Figures S4 and S5 (PC was used here for connectivity calculation).Other segmentation methods had similar results for the grid search on (W r , p KDE ) .Since future application of our method involves studying dynamic changes in the functional connectivity networks and using network features as inputs for neural decoding, the number of segments should not be too small.We also chose parameters so that the number of detected segments is similar to the number of sound stimuli, to have comparable p diff values.

Segmentation results
The EEG data of 28 participants were segmented and the difference ratio p diff values are shown in Fig. 9.For both Pearson correlation (Fig. 9a) and phase-locking value (Fig. 9b) connectivity matrices, graph distance measures DMDC, DMCC, DMNNC, and DMCCO performed better than the other segmentation methods in terms of p diff , meaning that intra-segment differences were smaller than inter-segment differences.The segmentation based on sound stimulation (SS) performed poorly, suggesting that the sound stimulus onset times did not mark boundaries in network structure change.This concurred with previous decoding analysis results on this data, which found that windows around the sound stimulation times were not the best choice of epochs for mental workload classification 51 .For the same length window segmentation, network structure between segments was not significantly different, showing that this naive approach is inappropriate.For GD, data from 4 participants' easy condition were not segmented, and data from 2 participants' difficult condition were not segmented (i.e.no boundaries were found).These unsegmented participants are not included in Fig. 9, but the numbers of segments are included in Supplementary Figure S6.For GD, the data from the participant with highest p diff in the difficult condition (value near 1.0) were divided into only two segments.The number of segments for each method is detailed in Supplementary Figure S6.From Supplementary Figure S6a and b, we can see that the number of segments found by KS varies greatly among participants, and the number of segments for other methods had smaller variance.Overall, DMDC, DMCC, DMNNC, and DMCCO outperformed (as in, inter-segment differences were larger than intra-segment differences) the methods based on matrix difference measures (GD and KS), as well as the two naïve methods for segmentation (SS and SL), on real data.These results support the higher effectiveness of the functional connectivity structure approach to segmentation when used on real EEG data.
In Figs. 10 and 11, we visualize the segmentation performance of the segmentation methods with graph distance measures.EEG data (28 × 2 data) and segmentation methods corresponding to the highest p diff are selected to visualize.The node importance over time is obtained using the sliding window method.We found that the fluctuations of EC and BC of each channel over time are relatively more disordered, which is not conducive to segmentation, the p diff value using DMEC and DMBC is low in Fig. 9.For other segmentation methods using graph distance measures, the boundary is close to the inflection point on the curve.

Discussion
In this study, we proposed an EEG segmentation method based on graph distance measures to solve the problem of how to segment EEG data according to changes in the functional connectivity network structure so that within each segment, differences in network structure remain small while between segments, differences are large.We designed a graph distance framework to provide an intuitive scalar measure for differences in network structure.We combine graph distance with a segmentation procedure that looks for abnormally large distances between a sliding window and a reference window.We evaluated our segmentation method on simulated EEG data under two network structure changes.Finally, we segmented real EEG data and quantified the connectivity changes within segments versus between segments.Although the closeness centrality graph distance measure (DMCC) performed well on both simulated data and actual data, its computational cost was relatively high.The main reason is that the calculation of shortest paths for closeness centrality has high computational cost.We used Dijkstra's algorithm 54 to calculate the shortest path between nodes.Considering that our graphs are dense, i.e. we need to traverse all edges, the time complexity of CC is O(n 3 ) , where n is the number of EEG channels.
Graph distance measures DMDC and DMCC have relatively good segmentation performance, as seen in Figs. 5, 6, and 9. DMBC performed well on Simulation 2 but poorly on Simulation 1. DMNCC performed poorly on both simulations but performed well on the real data.The six graph distance measures use node indices that ascribe importance to nodes in different ways, which leads to varying sensitivity to different network structures.In real applications, which measure is best suited may be difficult to predict a priori; therefore, we recommend the use of multiple graph distance measures and performing model selection.Generally, there are many indices for importance of nodes, such as Katz's centrality 55 , PageRank 56 , LeaderRank 57 , and semilocal centrality 58 .The choice of node importance index requires attention to interpretation and computational complexity.We can select the node importance index according to the network structure changes of interest.For example, closeness centrality describes the transmission efficiency of this node (brain region) to all other nodes.Using DMCC, at the time points of segmentation, there are large changes in the transmission efficiency of brain regions.There are some relationships between the six node importance indices.They characterize different aspects of network structure, and the relationship between them thus depends on the specific network characteristics or structure 59 .If there is prior knowledge about the form of the graph of the functional connectivity network, then we can say something more specific, for example, in scale-free networks, DMDC and DMBC are strongly correlated 60 .The flexibility of our approach allows different indices to be tested for an application.How indices relate to each other and the theory behind which index is best for segmentation in a given situation are important questions, which we do not try to answer here; but we have presented empirical comparisons in our simulations and analysis on real EEG data.
In our segmentation procedure, the length of the reference window W r and the probability p KDE in the detection of outliers have large impacts on the number of segments.The reference window determines the minimum length of detected segments.A shorter reference window can be used if one believes changes in EEG data happen frequently.The larger the probability p KDE , the less sensitive the outlier detection, and the number of detected segments will be smaller.When one believes the network structure has slow changes, we recommended a smaller p KDE .
The KS matrix measure represents previous work, the Source-Informed Segmentation method (SIS) 23 .Many steps of the presented sliding-window change-detection mechanism are also based on SIS's sliding-window algorithm.However, there are crucial differences in the overall framework.In the SIS method, the patterns of EEG voltage values of multiple channels (the EEG data matrix) are examined for changes, while in our framework, the differences in a chosen graph measure computed from functional connectivity networks (e.g.connectivity matrices) built from EEG data are examined.Our method also differs in the change detection mechanism: we identify a time window that has a distance value (or p-value) which is an outlier of the reference distribution of distances (p-values).In contrast, SIS 23 uses the p-values from K-S testing to directly decide significance.Our approach prevents situations where the largest distance (or smallest p-value) is not necessarily significantly higher (lower) than the other distances (p-values).
In Simulations 1 and 2, the GD method performed better than some graph distance measures.Especially in Simulation 2 at higher SNR, GD performed the best.Segmentation based on graph distance measures are relatively more sensitive to small network changes.In addition, the results of GD and DMEC were found to be coincident in the two simulations (Figs. 5, 6).Their results on real data are also similar (Fig. 9).This is likely because they both are calculated by eigenvectors, but we do not have an analytical proof of equality.
There are some limitations to our study.There exist other segmentation methods which we did not compare here.In this study, we only compare to GD and KS (as well as to stimulation-based segmentation and same-length window segmentation), because these also look at differences in functional connectivity, while other methods in the field optimize for other criteria.Each segmentation method has a different definition or goal and thus may be better suited for a particular analytical purpose.Another limitation is that our real data experiments were conducted using electrodes as network nodes, which means computed networks may include artifactual edges due to electrical conduction.However, this noiseness in the data is faced by all methods in our experiments.
We only consider single-layer undirected weighted networks.When there is a negative correlation, we use the absolute value for the weight.As a result, some information will be lost.In follow-up work, negative correlation can be considered and more complex networks can be constructed, such as directed networks or multi-layer networks.We only use one simulation method to make our simulated data, so performance under other simulation methods needs further verification.
We only tested our method with correlation and PLV as the functional connectivity metric, while many other metrics exist.Our framework is compatible with other metrics, but whether our method will work well with other metrics is an open question, and likely depends on the application and data characteristics.

Conclusion
In this work, we designed a graph distance framework and a segmentation procedure for the segmentation of EEG data.According to our results on simulated EEG data and real EEG data (with correlation or PLV as connectivity metric), the graph distance measures based on degree centrality and closeness centrality outperform matrixbased measures (Grassman distance and Kolmogorov-Smirnov).When also considering computational cost, the degree centrality approach presents a good overall choice.Performance of various node importance indices vary per dataset, suggesting that selection of node importance index is dependent on characteristics of data.
We provide a new perspective for segmentation, one based on network structure changes, and demonstrate its effectiveness with correlation and PLV connectivity, which suggests more segmentation methods based on network structure properties can be investigated.Our results suggest that when using an appropriate graph distance measure for the data, segmentation can be more sensitive and effective than matrix-based segmentation.

Figure 1 .
Figure 1.Example segmentation detected via our method (DMDC feature, edges from correlation).8-channel EEG data (lower panel) and functional connections of three segments (upper panels).Only edges above 0.2 are shown to improve visualization.Size and color of the nodes indicate nodes' degree centrality measure.Color of edges represents Pearson correlation coefficient(absolute value), which is unitless.

Figure 2 .
Figure 2. Flowchart of procedure for calculating graph distance measures (proposed method) and matrix measures (comparison method).(A) Generating connectivity matrices from two windows of EEG data.(B) Comparing connectivity matrices directly (comparison method).(C) Comparing networks using Euclidean distance on node importance index values (proposed method).

Figure 3 .
Figure 3. Flowchart of the segmentation algorithm.

Figure 4 .
Figure 4. Illustration depicts changes in community structure in the two simulations.(A) Conceptual diagram of community changes in Simulation 1, wherein a new community appears (as shown) and subsequently disappears.Left side shows weighted adjacency matrix; right side shows corresponding weighted network.Dashed lines indicate the regions where changes have occurred.Width and color of edges in the network are determined by weight values.Communities are distinguished by the color of nodes.(B) Conceptual diagram of changes in Simulation

Figure 9 .
Figure 9. Segmentation performance as measured with p diff on real EEG data (higher values are better). (a) connectivity calculated by Pearson correlation, (b) phase-locking value.SS denotes segments with boundaries at sound stimulation times.SL denotes segmentation into same-length time windows.

Figure 10 .
Figure 10.Visualization of results of segmentation based on different node importance indices on networks constructed by Pearson correlation.EEG data and segmentation results of the highest p diff are shown.Curves in the graph show node importance of all nodes, obtained in sliding windows (length 2 s, overlap 1 s).The data of the first 100 s are displayed here.Vertical dashed lines are segment boundaries.

Figure 11 .
Figure 11.Visualization of results of segmentation based on different node importance indices on networks constructed by PLV.EEG data and segmentation results of the highest p diff are shown.Curves in the graph show node importance of all nodes, obtained in sliding windows (length 2 s, overlap 1 s).The data of the first 100 s are displayed here.Vertical dashed lines are segment boundaries.

Table 2 .
Segmentation algorithm parameters.Initialize: length of the reference window W = W r .Number of possible boundaries k = 0 .List of distance values L = {} .Start of reference window t i is set to the boundary detected previously (or start of data).The reference window starts from time t i and goes to time t i + W r .The sliding window starts from time t i + W r − W v and goes to time t d Number of comparisons before starting outlier calculation ( W KDE d ≤ W d )

Table 3 .
SNR corresponding to signal strength k s values.