Computational Intelligence and Neuroscience

Volume 2018, Article ID 2406909, 15 pages

https://doi.org/10.1155/2018/2406909

## Assessment of Multivariate Neural Time Series by Phase Synchrony Clustering in a Time-Frequency-Topography Representation

Correspondence should be addressed to M. A. Porta-Garcia; moc.liamg@atropam

Received 25 December 2017; Accepted 30 January 2018; Published 21 March 2018

Academic Editor: Victor H. C. de Albuquerque

Copyright © 2018 M. A. Porta-Garcia 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.

#### 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 Time-Frequency-Topography 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 steady-state artifacts at 7.6 Hz. Further, contrast maps of Levenshtein Distance highlight synchrony differences between ERP and no-ERP 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 Time-Frequency (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 event-related 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* event-related 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 Cross-Coherence (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 Random-Matrix 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 equal-time 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 Time-Frequency-Topography (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 distribution-Rihaczek (RID-Rihaczek) 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 mean-field 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 one-to-one 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 Brain-Computer 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 across-trials 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 channel-wise. 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 no-ERP 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 discrete-time samples are weighted by the matrix , which is determined by the inverse-square 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://bnci-horizon-2020.eu/database/data-sets). The subjects were patients with amyotrophic lateral sclerosis and were naive to BCI training. The authors recorded eight EEG signals according to 10-10 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 band-pass-filtered 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 four-quadrant 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 so-called* 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