Recent Advances in Brain Signal Analysis: Methods and Applications 2018
View this Special IssueResearch Article  Open Access
Assessment of Multivariate Neural Time Series by Phase Synchrony Clustering in a TimeFrequencyTopography Representation
Abstract
Most EEG phase synchrony measures are of bivariate nature. Those that are multivariate focus on producing global indices of the synchronization state of the system. Thus, better descriptions of spatial and temporal local interactions are still in demand. A framework for characterization of phase synchrony relationships between multivariate neural time series is presented, applied either in a single epoch or over an intertrial assessment, relying on a proposed clustering algorithm, termed Multivariate Time Series Clustering by Phase Synchrony, which generates fuzzy clusters for each multivalued time sample and thereupon obtains hard clusters according to a circular variance threshold; such cluster modes are then depicted in TimeFrequencyTopography representations of synchrony state beyond mere global indices. EEG signals from P300 Speller sessions of four subjects were analyzed, obtaining useful insights of synchrony patterns related to the ERP and even revealing steadystate artifacts at 7.6 Hz. Further, contrast maps of Levenshtein Distance highlight synchrony differences between ERP and noERP epochs, mainly at delta and theta bands. The framework, which is not limited to one synchrony measure, allows observing dynamics of phase changes and interactions among channels and can be applied to analyze other cognitive states rather than ERP versus no ERP.
1. Introduction
There is a growing interest among the neuroscientific community to unravel the intricate neural mechanisms involved in the broad integration of different brain structures, which enable the emergence of cognitive processes. Several studies conducted with electroencephalography (EEG) and magnetoencephalography (MEG) have provided evidence that supports the idea of neural synchronization intrinsic to mental tasks, with the fluctuating disposition of communication channels in the nervous system, especially between active regions in the brain [1–5].
In this regard, phase locking analysis of neural oscillations and other different measures of synchronization has gained attention, as several methods have been developed to provide a quantitative view of synchronism in brain sources and their behavior, estimating phase synchrony (PS) from different perspectives, depending on the purpose of the study in question [6]. This same variety of methods and proposals causes lack of agreement in the terminology used to refer to all these measures. Roach and Mathalon have provided a wide review attempting to clarify this situation [7]. Thus, for the sake of following a standard of terms, descriptions of any PS measure will follow the referred publication.
In order to perform PS analysis, instantaneous phase information of EEG signals must be extracted. Most methods are based on wavelet analysis [6–10]. Another common technique besides wavelets for extracting instantaneous phase values from the analytical signal is the Hilbert transform. Analytic phase from wavelets or Hilbert transform has been shown to give almost same results as Short Time Fourier Transform adjusting the filter settings adequately [11, 12]. There are also other TimeFrequency (TF) decompositions used for obtaining phase information, such as Rihaczek distribution, Empirical Mode Decomposition, sinusoidal quadrature filters, and Matching Pursuit [13–16].
In general, for the study of PS, it can be said that there are two main approaches: phase locking and phase coherence. The former refers to the eventrelated phase locking across trials regarding an event’s onset over one electrode, that is, the Phase Locking Factor (PLF). If instantaneous phase angles between trials are closer to a uniform distribution over the unit circle, the PLF is close to zero; otherwise, it is close to one if instantaneous phase angles between trials are highly synchronized in the same direction over the unit circle. The latter approach, phase coherence, also called Phase Locking Value (PLV), or within the context, the eventrelated phase coherence across trials, evaluates consistency of phase differences between 2 electrodes across trials, also with values between 0 and 1. As the reader already noticed, each measure determines different types of PS; therefore, both measures can be complementary to each other [7].
Other types of measures, such as linear coherence or magnitude squared coherence, are not suitable to analyze PS; unlike PLF and PLV, both measures yield results weighted by magnitude, and the interpretation of these becomes unclear, since phase synchronization patterns and amplitude changes are not necessarily related to the same neural process [6–8]; Rosenblum demonstrated that PS of chaotic oscillators is possible, where bounded phase differences exist and variations of amplitude are chaotic and uncorrelated [17]. The Phase CrossCoherence (PCC) eliminates amplitude information and produces a function of phase differences averaged across trials [6].
All PS measures mentioned above focus on the evaluation of intertrial phase consistency over an individual EEG channel or phase differences between signals from two recoding sites, that is, providing only univariate or bivariate approaches. Nevertheless, the complete scenario involves a multichannel recording; thus a bivariate approach may not capture relevant information of all the dynamics and interactions of the full system [18, 19]. Thereupon, existing methods of multivariate synchronization analysis comprise even other metrics besides PS, based on different types of correlation measures. Correlation between probabilities of recurrence is used to measure PS, clearly distinguishing preseizure and seizure states of epileptic EEG [20, 21]. Based on RandomMatrix Theory (RMT), Osorio and Lai compute the average phase synchronization times (APSTs) among pairs of channels in order to construct a matrix, from which they use both the determinant and the eigenvalue spectra for assessing synchronization [22]. Li et al. presented another method based on RMT, using equaltime correlation instead of PS, and then the eigenvalue decomposition is used to calculate a global synchronization index that increases during epileptic seizures [23]. Mutlu et al. extend the concept of phase differences between two signals, mapping these differences onto an dimensional hyperspherical coordinate system; however, the authors later reported that Hyperspherical Phase Synchrony (HPS) is dependent on how the phase differences are sorted, which is corrected with another hyperdimensional coordinate system [19, 24]. Alba et al. proposed a visualization system with multitoposcopic graphs and TimeFrequencyTopography (TFT) maps for synchrony patterns, indicating increase, decrease, or an equal level of synchronization between pairs of electrodes with respect to a previous state, using different PS bivariate measures [15].
Some other approaches aim to improve the resolution of the TF decomposition used for extracting phase information. Aviyente and colleagues used a reduced interference distributionRihaczek (RIDRihaczek) for computing PLV [25]. Subsequently, the authors extend their method to quantify all possible pairwise comparisons and analyze those interactions between electrodes through a graph clustering algorithm, which allows overlapping clusters, and each electrode has a “participation score” that reflects their significance in the formation of a cluster [26]. Previous works also conceive the idea of clustering with degrees of membership. Allefeld and Kurths addressed the multivariate synchronization as a meanfield cluster of oscillators that participate in different degrees, that is, how close an oscillator phase is close to a reference phase, which is determined by the circular mean of all oscillator phases [27]. Nevertheless, the single cluster assumption dismisses other possible cluster formations. Later, the authors made a generalization of the cluster analysis to correct this issue based on eigenvalue decomposition of a matrix containing indices of bivariate synchronization strength, associating each eigenvalue greater than one to a cluster [28]; however, the onetoone correspondence between dominant eigenvectors and clusters is not always fulfilled [29].
Summarizing, multivariate methods help in perceiving overall synchronization patterns, providing a global index instead of matrices of bivariate comparisons [19]. Since many of these investigations focus on epilepsy studies, it makes sense to provide a general assessment of the synchronization state of the system with a crisp numerical value in order to distinguish seizure and preseizure conditions. Rather than a global index and aiming to characterize a broader variety of cognitive states, such as mental tasks for BrainComputer Interface (BCI), the framework proposed in this article points to observing the dynamics of phase changes along multivariate neural time series over the TF plane and projecting their interactions in TFT maps.
The proposed clustering algorithm, Multivariate Time Series Clustering by Phase Synchrony (mCPS), establishes local relations by means of clusters of highly synchronized signals in each sample time, allowing exploring these phase associations through all samples searching for patterns of cluster formations. Additionally, our proposal also addresses an acrosstrials perspective. Thus, it can be said that the PS measure used in this work is more related to PLF (circular variance) rather than to phase coherence (consistency of phase differences), applied channelwise. Haig et al. proposed a similar conception of PS, which lacks an automatized selection of synchronized signals via clustering [30].
Beyond yielding a PS measurement and a TFT portrayal, the framework also provides contrast maps of Levenshtein Distance (LD) as a metric for visual analysis and comparison of differences in PS patterns between different conditions (in this case, ERP and noERP epochs), as well as TF images of channels, highlighting which clusters of PS can be related to the changes of power due to the ERP. While some of the methods mentioned before use clustering analysis, like [26], most of them are fuzzy clusters in short time windows and without topographic representation. The way mCPS is conceived requires hard clustering, as it will be further detailed.
2. Materials and Methods
2.1. Simulated EEG and Experimental Data
Several experiments were carried out with both synthetic and real EEG signals (sEEG and rEEG, resp.) in order to determine the extent to which our framework is capable of retrieving reliable and useful information (presented as clusters of electrodes) that allows establishing relationships between highly synchronized EEG channels and the brain activity of interest through time samples and over different bandwidths. The sEEG was built based on a linear mixing model of independent sources , with a sampling frequency of 256 Hz, resulting in observed signals . Contributions of every through the discretetime samples are weighted by the matrix , which is determined by the inversesquare law of distances between and locations:
Spatial location of each electrode corresponds to the basic 10–20 international system [32] over a unit sphere. The volume conduction of the EEG model was assumed to be homogeneous and isotropic. The complete sEEG record is constructed with 30 epochs of 3 seconds, each of them containing a simulated Visual Evoked Potential (VEP) centered at 1.5 s from the epoch onset (peak amplitude at 1500 ms and constant across trials). Equation (2) describes the construction of the VEP:where Hz, , and . Besides the VEP, sources comprise three different types of noise components: (a) harmonics, which vary in amplitude, frequency of the sinusoidal oscillations, and initial phase and (b) white Gaussian and (c) colored Gaussian noise. Localization of within the brain area of the model can be either a fixed position or a linear displacement or with rotational motion.
In order to assess the framework with rEEG, four subjects (S2, S5, S6, and S7) were selected from a record of P300 evoked potentials [33] using the P300 Speller proposed by Farwell and Donchin [34] (available at http://bncihorizon2020.eu/database/datasets). The subjects were patients with amyotrophic lateral sclerosis and were naive to BCI training. The authors recorded eight EEG signals according to 1010 standard (Fz, Cz, Pz, Oz, P3, P4, PO7, and PO8) using active electrodes, referenced to the right earlobe and grounded to the left mastoid. EEG signal was digitized at 256 Hz and bandpassfiltered between 0.1 and 30 Hz. Subjects were required to spell seven predefined words of five characters each by controlling the P300 matrix speller. It should be mentioned that no extra preprocessing was performed over the data. The first three runs (15 trials in total) are described as “calibration runs” and runs 4–7 are the “testing runs” where participants were provided with feedback.
2.2. Clustering EEG Channels according to Circular Variance
2.2.1. Extraction of Phase Information
Given (for sEEG, and, for rEEG, ) signals, a TF decomposition is performed over the continuous EEG with predefined bandwidths at center frequencies:where , Hz, and for both sEEG and rEEG. Such decomposition is carried out with a Continuous Wavelet Transform (CWT) at peak frequencies from (3) with complex Morlet wavelets:where is the standard deviation of the Gaussian function used to make each and is the number of wavelet cycles (in this case, ). Then, the instantaneous phase is obtained from (5), using implementation of the fourquadrant inverse tangent:
2.2.2. Multivariate Time Series Clustering by Phase Synchrony (mCPS)
Algorithm 1 explains how mCPS works, which is based on directional statistics to measure the degree of phase locking and formation of clusters. The circular spread in angular data can be computed with the magnitude of the socalled mean resultant vector [35]. Directional data (in this case, of the signals) can be observed as points over the unit circle. Then, the Cartesian coordinates of the center of mass can be expressed as , where

Therefore, . Magnitude of is close to 1 when EEG channels are highly phaselocked; it is close to zero otherwise. PortaGarcia et al. presented an example using magnitude changes of vector over time in a determined group of EEG channels comparing two different conditions [36]. The functioning of mCPS over EEG channels according to circular variance is as follows.
Once is retrieved for the entire EEG, the procedure CreateFuzzyClusters generates fuzzy clusters of electrodes for each time sample and for each center frequency . The threshold defines whether or not an electrode is assigned to a determined , and as fuzzy clusters consider intersections of cluster elements, the main task of the procedure ConvertToHardClusters is to obtain hard clusters by preserving clusters with higher value of and proceed to eliminate intersections iteratively of the remaining in such a way that . Therefore, the result of mCPS is a new matrix cEEG, containing the cluster labels to which each EEG channel belongs in each time sample .
2.2.3. Cluster Labeling
Every run of mCPS is bounded for each time sample , and an arbitrary numeric label is assigned to each cluster. Then, an example of generated clusters could be for and for . In this case, numeric labels 1 and 2 do not provide any useful information of cluster content. In order to establish a meaningful relationship that reflects that and are actually the same cluster, a labeling system was developed based on hexadecimal words that encode which electrode belongs to the cluster and then assign a specific color in a onetoone relationship to represent clusters in a TFT map, which will be described further. In Figure 1, it can be observed that each hexadecimal digit corresponds to binary bits of electrode quartets, where digit 1 means that the electrode is assigned to a determined cluster if magnitude of is greater than threshold . Therefore, a hexadecimal word of two digits encodes the cluster membership for eight EEG channels. As a consequence of this encoding system, a slightly different hue of color label should depict similarity between clusters, for example, a blue cluster containing electrodes P3, Pz, and P4 and a lighter blue cluster that only contains P3 and Pz. Hence, the matrix cEEG now has as elements the cluster labels of hexadecimal words.
2.2.4. Construction of TimeFrequencyTopography (TFT) Maps
To be able to condense the large amount of information obtained from mCPS and make it suitable for visual analysis, we used TFT maps for topographic representation of all yielded in previous steps. Some previous schemes of TimeFrequencyTopographic visualization can be found in literature [15, 37]. Then, the cEEG section that corresponds to the rEEG segment to be analyzed is windowed, displaying scalp maps with cluster modes of the cEEG windows of size , which is specified in number of samples (Figure 2). The cluster modes for each channel are assigned only if the mode frequency is greater than threshold . For both rEEG and sEEG, ; this way, bimodal or multimodal results are avoided.
With regard to the rEEG, it should be mentioned that as the selected runs for analysis with our method were clustered separately, the color labels in a TFT map of ERP condition are the same as a TFT map of no ERP only if it is the same subject and same run; otherwise, this condition may not be satisfied, except for two cases: the color map is bounded between specific RGB values between dark blue and bright yellow, which corresponds with cluster with hexadecimal label “01” (only channel P8 is assigned) and cluster “FF” (all channels are assigned), respectively. Intermediate variations of label color depend on the amount of generated clusters along time.
2.2.5. Intertrial TFT (iTFT) Maps
An iTFT depicts modes across epochs. It can be seen as a TFT map containing intertrial cluster modes (ITCM) instead of computing cluster modes over a cEEG segment directly (Figure 3). Regarding the rEEG, for each run of the experimental protocol, the instantaneous phase is computed over the complete run and the clustering is performed before epoch segmentation. After these steps, ERP and noERP epochs are taken separately and their ITCM is computed in such a way that the most representative cluster formations over the ERP and noERP epochs are retrieved. For the rEEG case, the resultant iTFT map illustrates the most prevalent phase clustering patterns over 1000 ms (duration of trials) with a time window of size (62.5 ms).
2.2.6. Levenshtein Distance (LD) and Complementary TF Maps
LD is included to sense relevant differences between ERP and noERP epochs. This measure can be defined as the minimum cost of transforming one string into another through a sequence of operations [38]:where and are strings constructed with characters of the same alphabet and represents the set of edit operations to make , weighted by function . With being a simple edit operation and being the null string, there are three types of transformations: insertions , substitutions , and deletions . Adapted to our case, , , and and are binary cluster labels of same length; thus the only operation to perform is substitutions of characters. Since clusters labels encipher the membership of electrodes, the maximum LD should be equal to 8 for the extreme case of “00000000” (which means that no cluster mode was assigned to any channel due to threshold ) and a cluster mode with all 8 electrodes within .
Furthermore, additional TF maps are generated from the CWT of each channel, which coupled with LD measures, and they help to observe findings in the mCPS information that could be associated with the changes of power due to the ERP over the time series. The LD distances are depicted in TimeFrequencyLevenshtein (TFL) maps.
2.2.7. Framework Pipeline
The complete framework pipeline is shown in Figure 4. Once the extraction of phase information of EEG in block and mCPS is performed in block , EEG clusters (cEEG) are labeled in block and then segmented according to the acquisition protocol. For this particular case, condition 1 and condition 2 in Figure 4 correspond to ERP and noERP epochs, respectively. Important to notice, segmentation of cEEG occurs after the hexadecimal labeling (block c) in order to allow direct comparison between conditions in the iTFT maps, ensuring a onetoone correspondence among color labels in the topographic scalp layouts of clusters and hexadecimal labels. Finally, the TFL maps (block ) highlight dissimilarities over time and frequency of the mCPS outputs for ERP and no ERP.
3. Results
Figure 5 summarizes the most remarkable outcome of the experiments with sEEG. Figure 5(a) shows the grand average of each channel, and Figure 5(b) displays the corresponding spectra of all channels as well as the scalp distribution of power at center frequency Hz. Figure 5(c) shows a single scalp map extracted from the correspondent TFT maps after applying mCPS over a single trial of sEEG, positioned at 1500 ms (which is where the peak amplitude of the VEP is found) and centered at , with a signaltonoise ratio (SNR) of 0.328 dB. Figure 5(d) also shows a single scalp map, at same latency and center frequency , coming from a TFT map generated after applying mCPS over the grand average of the 30 epochs, with SNR = 3.16 dB. By visual inspection, it can be observed in Figure 5(d) that electrodes in blue cluster correspond to those in Figure 5(a), where the VEP is more evident (marked with red circles); it also largely coincides with the scalp areas with highest power at (Figure 5(b)). Remarkable to say, despite the lower SNR in a single trial compared to scalp map of Figure 5(d), mCPS is able to retrieve some of the electrodes within the blue cluster (Figure 5(c)).
(a)
(b)
(c)
(d)
With respect to rEEG, the main attention was on the intertrial analysis searching for differences between ERP and noERP conditions, using iTFT maps. Different values of threshold were tested between 0.90 and 0.99 for cluster mode assignments, while was fixed at 50% and samples. In relation to the data, from the seven runs of each subject, only the testing runs (4–7) were processed with our framework, each of them individually. For reasons of space, only some relevant portions of maps per subject are presented in figures: run 4 for S2 and for S5 ( and , resp.), run 7 for S6, and run 6 for S7 (both with ). Figure 6 shows grand averages of all channels for these runs for each subject, contrasting ERP condition (blue) versus noERP condition (red). Respecting the TFL and TF maps, only the most illustrating channel is depicted. For the full maps of the runs mentioned before, please refer to http://itzamna.uam.mx/lini/mcps.html.
Results of run 4 for S2 are displayed in Figure 7. In the ERP iTFT map (Figure 7(a)), formations of cluster modes with label “FF” (bright yellow) containing P3P4PO7OzPO8 can be observed from 312.5 ms to 750 ms at 2 Hz. The same situation occurs at 2.5 Hz with P3PO7P4. No characteristic cluster formation is shown in the noERP iTFT map (Figure 7(b)). Noteworthy, run 5 portrayed similar conditions compared to run 4, except that relevant cluster formations were found in bins centered at 1.3, 1.6, and 2 Hz. As for runs 6 and 7, neither ERP nor noERP iTFT maps of S2 revealed any characteristic cluster formation. In Figure 7(d), the TFL map for P4 is displayed. Important to highlight, this map depicts yellow areas that coincide (at least visually) with the concentration of power of the P300 wave (Figure 7(e)), particularly for P3, P4, PO7, and Oz (TFL maps for P3, PO7, and Oz can be observed in the complete study). It is also coincidental with the cluster formations described previously in the ERP iTFT map and with the P300 power time course, around 312 ms and 750 ms approximately (Figure 7(e)), which is not the case if such cluster arrangements are compared with noERP TF maps (Figure 7(f)).
(a)
(b)
(c)
(d)
(e)
(f)
3.1. SteadyState Visual Evoked Potential (SSVEP) Artifact
As depicted in Figure 7(c), it can be observed that cluster “FF” contains all EEG channels at 7.6 Hz over the entire row. This is highly likely to be related to an SSVEP artifact derived from a fixed value of the interstimulus duration (125 ms). This pattern appears in all subjects, with some minor variations of . For example, for S2, this fact can be related to the concentration of power around 7.6 Hz in the entire epoch in all TF maps of each channel (Figures 7(e) and 7(f)). This can be verified with almost all TF maps presented for both ERP and no ERP for all subjects.
Regarding S5, observations within runs are very similar (). For ERP condition (Figure 8(a)), cluster formations of “FF” with parietal channels and Cz clearly coincide with yellow areas of TFL map of P3 (Figure 8(c)) and power concentration of P300 in the TF map (Figure 8(d)). This can be observed over the ERP iTFT cluster formations in bins centered at 1.6, 2, 2.5, and 3.1 Hz (Figure 8(a)). iTFT maps of no ERP (Figure 8(b)) did not show any relevant cluster formation.
(a)
(b)
(c)
(d)
(e)
For S6 (), Figures 9(a) and 9(b) illustrate a section of the correspondent iTFT map of ERP and noERP conditions, respectively. The “FF” cluster formations can be observed in the 3.1 Hz bin, which takes place at different time windows. There are no relevant cluster formations over noERP map at the same times. The TFL map confirms these differences with the yellow areas for Pz (Figure 9(c)). In this case, the relationship with the power in TF map of P300 wave (Figures 9(d) and 9(e)) is not so evident.
(a)
(b)
(c)
(d)
(e)
Concerning S7, in run 6 with a threshold of , parietal electrodes stand out again portraying diverse “FF” cluster arrangements over the scalp, mainly at 2, 2.5, and 3.1 Hz (Figure 10(a)), concurring with yellow areas in corresponding TFL map of P3 (Figure 10(c)) and with the power of P300 wave in the TF map (Figure 10(d)). In run 7, the appearance of other cluster formations besides “FF” (perceived in other runs and subjects) was noticeable, with parietal electrodes between 375 and 625 ms.
(a)
(b)
(c)
(d)
(e)
4. Discussion
The findings over the TFT maps of sEEG served as a starting point for leading the research to the analysis with real data, as coincidences of the generated cluster in the single trial and the one over the grand average reflected the ability of mCPS to retrieve the PS information of interest. For rEEG, the iTFT maps exposed several differences between ERP and no ERP, and maybe the most notable and consistent feature was the arrangements of clusters labeled “FF” systematically appearing in ERP maps (Figures 7(a), 8(a), 9(a), and 10(a)) derived from the ITCM computation, contrasting with the absence of such patterns in the noERP maps (Figures 7(b), 8(b), 9(b), and 10(b)). This fact is evidenced with the TFL maps (Figures 7(d), 8(c), 9(c), and 10(c)), highlighting the areas of the TF plane with perceptible differences among ERP and no ERP. Such differences can be noted by contrasting Figures 7(e), 7(f), 8(d), 8(e), 9(d), 9(e), 10(d), and 10(e), respectively. Moreover, most of the “FF” appearances can be related (at least by visual inspection) to the P300 wave, given the times and bandwidths where these clusters appear, as most of them were localized within delta and theta ranges, which is consistent with frequency content of a P300 ERP [39–41]. The frequency content of noERP epochs observed in Figures 7(f), 8(e), 9(e), and 10(e) could hardly be explained by any neurophysiological event of relevance, but rather it could be due to subharmonics of the SSVEP artifact, as the power concentration can be perceived as extended “lines” throughout the time series.
There were cases (like S6) where analysis with TFL and TF maps did not yield any clear distinction between ERP and no ERP, like run 6, where cluster formations were sporadic and intermittent, making it difficult to establish a relationship with the P300 wave. Noteworthy, samples of ERP and noERP epochs are highly unbalanced (each run per subject contains 100 ERP epochs and 500 noERP epochs), which reinforces our results distinguishing these conditions, considering the fact that we are using mode as statistical measure, and despite a greater amount of samples of noERP epochs, no relevant cluster modes formations were detected.
Another important aspect is related to frequency locking and tracking of frequency flows [42, 43]. A limitation in some methods relying on a narrow band TF decomposition, such as the frequency bins generated with wavelets, is the assumption of frequency stationarity of PS, hiding or masking periods of continuous PS with transient variability of synchronization frequency through time. However, the TFT maps can capture this frequency flow of PS, as it can be observed how the “FF” cluster patterns appear over different low frequency narrow bands, such as = 1.6, 2, 2.5, and 3.1 Hz for S5 (Figure 8(a)) and = 2, 2.5, and 3.1 Hz for S7 (Figure 10(a)).
With respect to hyperparameters, further analysis should be made varying threshold values ( and ) in order to evaluate the produced effect in cluster generation and visualization. As mentioned before, several tests were made with different values of , yet the results shown in this work are only for one per subject, which was heuristically selected by identifying the TFL maps that yield a better differentiation of ERP and noERP conditions. The method is highly sensitive to variations, and future work can be directed to automate selection of optimal values for .
Even though in these results our framework serves in identifying PS dynamics related to the neural activity of interest organized and structured in clusters of EEG channels, there is still a lot of room for improvement. At this point, our method describes nearzero phase lag relationships between EEG channels ( in most of the cases). By definition, volume conduction requires zero phase lag, but a phase difference close to zero is not necessarily due to volume conduction, as this kind of phase associations can be found widespread over the cerebral cortex because of corticothalamic projections [44]. There are some measures such as Phase Lag Index (PLI) [45] or imaginary coherence [46] which deal with volume conduction by discarding zero phase lags, but at the same time these approaches are insensitive to true nearzero phase lag interactions [47].
On the other hand, volume conduction can be addressed by measuring phase reset, which can be detected when a phase shift takes place between two phaselocked signals [48]. This idea can be extended in our framework, trying to find phase resets between EEG channels. Adding other phase differences or phaselocking measures could retrieve different clustering patterns, which along with our already implemented mCPS measure and detection of phase resets could deliver complementary and relevant information.
5. Conclusions
Our framework provides a feasible way to address both single and intertrial PS analysis of multivariate neural time series, characterizing the PS variability through time. The majority of PS measures so far suggested in literature such as PLV or PCC are calculated between two signals [6–8] or provide only a global index of synchronization in the case of multivariate measures [19, 20, 23, 27]. Our framework is an alternative for studying the behavior of phase synchronization between all EEG channels at once in a given time window within different bandwidths of interest. Noticeable to say, the framework is not limited to any particular phase extraction technique (further discussion about the selection of these techniques is beyond the scope of this article) and can also easily be adapted to other PS measures like phase coherence, obtaining clusters of phase differences consistency from mCPS. It remains to assess and compare the proposed algorithm to other clustering algorithms in terms of efficiency and computational complexity.
The insight given by the iTFT maps provides a qualitative measure of intertrial cluster consistency, which when combined with the TFL and TF maps becomes helpful to assess which clusters patterns are related to a specific mental task. It should be mentioned that some yellow areas depicted in TFL maps that do not match with the power increase of the P300 wave shown in TF maps could be due to artifacts artificially derived from LD computation or due to other relevant neural information not related to ERP. Further analysis should be made regarding this issue.
Although in this first approach mCPS was applied over synthetic signals and P300 wave data with relatively few electrodes, the aim of this work was merely to illustrate the framework pipeline and how it describes PS patterns. As mentioned before, our work attempts to encompass a broader variety of cognitive states. For example, in the context of BCI, our framework might be useful for the characterization of mental tasks suitable for endogenous BCI paradigms with no external stimuli in the system. Then, feature extraction could be performed from mCPS outcome for asynchronous (selfpaced) BCI classification, distinguishing idle state from a specific mental task. Additionally, when exploring higher density EEG (64 channels or more), this framework could be used as a channel optimization tool finding the clusters of electrodes that contribute the most to characterization of a mental state.
Electrical signals from brain sources are volume conducted through nervous tissue, cerebrospinal fluid, skull, and scalp. Hence, an underlying issue in EEG recordings regards the single source contamination of multiple sensors via volume conduction. The EEG recorded over the scalp does not necessarily capture the direct activity underneath the electrode but a weighted mixture of different sources (neural or artifact). Then, distinction between volume conduction and true synchrony remains an open issue. Some authors have reported that methods for improving spatial resolution of EEG, such as scalp current density profiles (SCD), seem convenient as preprocessing steps before the estimation of PS [7, 8]. For future work, it should be interesting to study the effects of rereferencing. Again, in the BCI field, it could be assessed if rereferencing enhances performance using phase clusters as features for classification, bearing in mind the fact that the original phase delays may be distorted. It should be pointed out that no additional preprocessing was made, preserving the data as raw as possible. Further approaches for addressing volume conduction should be considered in forthcoming research.
Finally, to summarize the contributions, the proposed framework incorporates several features useful for PS analysis, such as iTFT and TFL maps, taking into account some aspects like frequency nonstationarity and flexibility of use of other synchronization measures besides PLF. The LD is applied as a metric for better distinction of differences between conditions, highlighting synchrony differences between ERP and noERP epochs, mainly at delta and theta bands. Additional information like the steadystate artifacts at 7.6 Hz is also retrieved and depicted in iTFT maps. Taking EEG as the view port of cortical activity, our framework provides a new insight into terms of largescale integration of emerging synchrony patterns of instantaneous phase during cognitive tasks, depicted in phaserelated cluster arrangements over the time series of EEG signals.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
M. A. PortaGarcia, M.S., gratefully acknowledges Consejo Nacional de Ciencia y Tecnología (CONACYT) for Ph.D. scholarship 271659.
References
 W. H. Miltner, C. Braun, M. Arnold, H. Witte, and E. Taub, “Coherence of gammaband EEG activity as a basis for associative learning,” Nature, vol. 397, no. 6718, pp. 434–436, 1999. View at: Publisher Site  Google Scholar
 J. Bhattacharya, H. Petsche, and E. Pereda, “Interdependencies in the spontaneous EEG while listening to music,” International Journal of Psychophysiology, vol. 42, no. 3, pp. 287–301, 2001. View at: Publisher Site  Google Scholar
 M. Popescu, A. Otsuka, and A. A. Ioannides, “Dynamics of brain activity in motor and frontal cortical areas during music listening: A magnetoencephalographic study,” NeuroImage, vol. 21, no. 4, pp. 1622–1638, 2004. View at: Publisher Site  Google Scholar
 J. Bhattacharya and H. Petsche, “Phase synchrony analysis of EEG during music perception reveals changes in functional connectivity due to musical expertise,” Signal Processing, vol. 85, no. 11, pp. 2161–2177, 2005. View at: Publisher Site  Google Scholar
 J. F. Cavanagh, M. X. Cohen, and J. J. B. Allen, “Prelude to and resolution of an error: EEG phase synchrony reveals cognitive control dynamics during action monitoring,” The Journal of Neuroscience, vol. 29, no. 1, pp. 98–105, 2009. View at: Publisher Site  Google Scholar
 M. S. Doesburg and M. L. Ward, “Synchronization between sources: emerging methods for understanding largescale functional networks in the human brain,” in In Coordinated Activity in the Brain, pp. 25–42, Springer, 2009. View at: Google Scholar
 B. J. Roach and D. H. Mathalon, “EventRelated EEG TimeFrequency Analysis: An Overview of Measures and An Analysis of Early Gamma Band Phase Locking in Schizophrenia,” Schizophrenia Bulletin, vol. 34, no. 5, pp. 907–926, 2008. View at: Publisher Site  Google Scholar
 J.P. Lachaux, E. Rodriguez, J. Martinerie, and F. J. Varela, “Measuring phase synchrony in brain signals,” Human Brain Mapping, vol. 8, no. 4, pp. 194–208, 1999. View at: Publisher Site  Google Scholar
 J.P. Lachaux, E. Rodriguez, M. Le Van Quyen, A. Lutz, J. Martinerie, and F. J. Varela, “Studying singletrials of phase synchronous activity in the brain,” International Journal of Bifurcation and Chaos, vol. 10, no. 10, pp. 2429–2439, 2000. View at: Google Scholar
 S. G. Roux, T. Cenier, S. Garcia, P. Litaudon, and N. Buonviso, “A waveletbased method for local phase extraction from a multifrequency oscillatory signal,” Journal of Neuroscience Methods, vol. 160, no. 1, pp. 135–143, 2007. View at: Publisher Site  Google Scholar
 M. le van Quyen, J. Foucher, J.P. Lachaux et al., “Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony,” Journal of Neuroscience Methods, vol. 111, no. 2, pp. 83–98, 2001. View at: Publisher Site  Google Scholar
 A. Bruns, “Fourier, Hilbert and waveletbased signal analysis: are they really different approaches?” Journal of Neuroscience Methods, vol. 137, no. 2, pp. 321–332, 2004. View at: Publisher Site  Google Scholar
 A. W. Rihaczek, “Signal Energy Distribution in Time and Frequency,” IEEE Transactions on Information Theory, vol. 14, no. 3, pp. 369–374, 1968. View at: Publisher Site  Google Scholar
 N. E. Huang, Z. Shen, S. R. Long et al., “The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis,” Proceedings A, vol. 454, no. 1971, pp. 903–995, 1998. View at: Publisher Site  Google Scholar  MathSciNet
 A. Alba, J. L. Marroquin, J. Peña, T. Harmony, and B. GonzalezFrankenberger, “Exploration of eventinduced EEG phase synchronization patterns in cognitive tasks using a timefrequencytopography visualization system,” Journal of Neuroscience Methods, vol. 161, no. 1, pp. 166–182, 2007. View at: Publisher Site  Google Scholar
 S. G. Mallat and Z. Zhang, “Matching Pursuits With TimeFrequency Dictionaries,” IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397–3415, 1993. View at: Publisher Site  Google Scholar
 M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Physical Review Letters, vol. 76, no. 11, pp. 1804–1807, 1996. View at: Publisher Site  Google Scholar
 S. Bialonski and K. Lehnertz, “Identifying phase synchronization clusters in spatially extended dynamical systems,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 74, no. 5, Article ID 051909, 2006. View at: Publisher Site  Google Scholar
 M. AlKhassaweneh, M. VillafaneDelgado, A. Y. Mutlu, and S. Aviyente, “A measure of multivariate phase synchrony using hyperdimensional geometry,” IEEE Transactions on Signal Processing, vol. 64, no. 11, pp. 2774–2787, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 D. Rangaprakash, “Connectivity analysis of multichannel EEG signals using recurrence based phase synchronization technique,” Computers in Biology and Medicine, vol. 46, no. 1, pp. 11–21, 2014. View at: Publisher Site  Google Scholar
 D. Rangaprakash and N. Pradhan, “Study of phase synchronization in multichannel seizure EEG using nonlinear recurrence measure,” Biomedical Signal Processing and Control, vol. 11, no. 1, pp. 114–122, 2014. View at: Publisher Site  Google Scholar
 I. Osorio and Y.C. Lai, “A phasesynchronization and randommatrix based approach to multichannel timeseries analysis with application to epilepsy,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 21, no. 3, 033108, 11 pages, 2011. View at: Google Scholar  MathSciNet
 X. Li, D. Cui, P. Jiruska, J. E. Fox, X. Yao, and J. G. R. Jefferys, “Synchronization measurement of multiple neuronal populations,” Journal of Neurophysiology, vol. 98, no. 6, pp. 3341–3348, 2007. View at: Publisher Site  Google Scholar
 A. Y. Mutlu and S. Aviyente, “Hyperspherical phase synchrony for quantifying multivariate phase synchronization,” in Proceedings of the 2012 IEEE Statistical Signal Processing Workshop, SSP 2012, pp. 888–891, USA, August 2012. View at: Publisher Site  Google Scholar
 S. Aviyente, E. M. Bernat, W. S. Evans, and S. R. Sponheim, “A phase synchrony measure for quantifying dynamic functional integration in the brain,” Human Brain Mapping, vol. 32, no. 1, pp. 80–93, 2010. View at: Publisher Site  Google Scholar
 M. E. Bolanos, S. Aviyente, and E. M. Bernat, “Identifying functional clusters in the brain using phase synchrony,” in Proceedings of the 2010 IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP 2010, pp. 5446–5449, March 2010. View at: Publisher Site  Google Scholar
 C. Allefeld and J. Kurths, “An approach to multivariate phase synchronization analysis and its application to eventrelated potentials,” International Journal of Bifurcation and Chaos, vol. 14, no. 2, pp. 417–426, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 C. Allefeld, M. Müller, and J. Kurths, “Eigenvalue decomposition as a generalized synchronization cluster analysis,” International Journal of Bifurcation and Chaos, vol. 17, no. 10, pp. 3493–3497, 2007. View at: Publisher Site  Google Scholar
 C. Allefeld and S. Bialonski, “Detecting synchronization clusters in multivariate time series via coarsegraining of Markov chains,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 76, no. 6, Article ID 066207, 2007. View at: Publisher Site  Google Scholar
 A. R. Haig, E. Gordon, J. J. Wright, R. A. Meares, and H. Bahramali, “Synchronous cortical gammaband activity in taskrelevant cognition,” NeuroReport, vol. 11, no. 4, pp. 669–675, 2000. View at: Publisher Site  Google Scholar
 A. Delorme and S. Makeig, “EEGLAB: an open source toolbox for analysis of singletrial EEG dynamics including independent component analysis,” Journal of Neuroscience Methods, vol. 134, no. 1, pp. 9–21, 2004. View at: Publisher Site  Google Scholar
 Trans Cranial Technologies, “Trans Cranial Technologies, 2410 Fortis Tower,” in 10/20 System Positioning Manual, pp. 77–79, Hong Kong, China, 2012. View at: Google Scholar
 A. Riccio, L. Simione, F. Schettini et al., “Attention and P300based BCI performance in people with amyotrophic lateral sclerosis,” Frontiers in Human Neuroscience, vol. 7, article 732, 2013. View at: Publisher Site  Google Scholar
 L. A. Farwell and E. Donchin, “Talking off the top of your head: Toward a mental prosthesis utilizing eventrelated brain potentials,” Electroencephalography and Clinical Neurophysiology, vol. 70, no. 6, pp. 510–523, 1988. View at: Publisher Site  Google Scholar
 K. V. Mardia and P. E. Jupp, Directional Statistics, John Wiley & Sons, New York, NY, USA, 1999. View at: Publisher Site  MathSciNet
 M. A. PortaGarcia, F. J. GarciduenasVargas, O. YanezSuarez, and R. ValdesCristerna, “EEG clustering based on phase synchrony for selfpaced BCI development,” in Proceedings of the In Proceedings of the 6th International BrainComputer Interface Meeting, 2016. View at: Google Scholar
 T. Koenig, F. MartiLopez, and P. ValdesSosa, “Topographic timefrequency decomposition of the EEG,” NeuroImage, vol. 14, no. 2, pp. 383–390, 2001. View at: Publisher Site  Google Scholar
 L. Yujian and L. Bo, “A normalized Levenshtein distance metric,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 29, no. 6, pp. 1091–1095, 2007. View at: Publisher Site  Google Scholar
 C. BasarEroglu, E. Basar, T. Demiralp, and M. Schürmann, “P300response: possible psychophysiological correlates in delta and theta frequency channels. A review,” International Journal of Psychophysiology, vol. 13, no. 2, pp. 161–179, 1992. View at: Publisher Site  Google Scholar
 V. Kolev, T. Demiralp, J. Yordanova, A. Ademoglu, and Ü. IsogluAlkaç, “Timefrequency analysis reveals multiple functional components during oddball P300,” NeuroReport, vol. 8, no. 8, pp. 2061–2065, 1997. View at: Publisher Site  Google Scholar
 T. Demiralp, A. Ademoglu, M. Schürmann, C. BasarEroglu, and E. Basar, “Detection of P300 waves in single trials by the wavelet transform (WT),” Brain and Language, vol. 66, no. 1, pp. 108–128, 1999. View at: Publisher Site  Google Scholar
 P. Tass, M. G. Rosenblum, J. Weule et al., “Detection of n:m phase locking from noisy data: application to magnetoencephalography,” Physical Review Letters, vol. 81, no. 15, pp. 3291–3294, 1998. View at: Publisher Site  Google Scholar
 D. Rudrauf, A. Douiri, C. Kovach et al., “Frequency flows and the timefrequency dynamics of multivariate phase synchronization in brain signals,” NeuroImage, vol. 31, no. 1, pp. 209–227, 2006. View at: Publisher Site  Google Scholar
 M. Steriade, “Grouping of brain rhythms in corticothalamic systems,” Neuroscience, vol. 137, no. 4, pp. 1087–1106, 2006. View at: Publisher Site  Google Scholar
 C. J. Stam, G. Nolte, and A. Daffertshofer, “Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources,” Human Brain Mapping, vol. 28, no. 11, pp. 1178–1193, 2007. View at: Publisher Site  Google Scholar
 G. Nolte, O. Bai, L. Wheaton, Z. Mari, S. Vorbach, and M. Hallett, “Identifying true brain interaction from EEG data using the imaginary part of coherency,” Clinical Neurophysiology, vol. 115, no. 10, pp. 2292–2307, 2004. View at: Publisher Site  Google Scholar
 J. M. Palva, “Phaselocking methods,” Encyclopedia of Computational Neuroscience, pp. 2356–2359, 2015. View at: Google Scholar
 R. W. Thatcher, “Coherence, phase differences, phase shift, and phase lock in EEG/ERP analyses,” Developmental Neuropsychology, vol. 37, no. 6, pp. 476–496, 2012. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 M. A. PortaGarcia et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.