{'fraction_non_alphanumeric': 0.029275808936825885, 'fraction_numerical': 0.011354746947967287, 'mean_word_length': 4.904268719384185, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identifying phase synchronization clusters in spatially extended dynamical systems

We investigate two recently proposed multivariate time series analysis techniques that aim at detecting phase synchronization clusters in spatially extended, nonstationary systems with regard to field applications. The starting point of both techniques is a matrix whose entries are the mean phase coherence values measured between pairs of time series. The first method is a mean field approach which allows to define the strength of participation of a subsystem in a single synchronization cluster. The second method is based on an eigenvalue decomposition from which a participation index is derived that characterizes the degree of involvement of a subsystem within multiple synchronization clusters. Simulating multiple clusters within a lattice of coupled Lorenz oscillators we explore the limitations and pitfalls of both methods and demonstrate (a) that the mean field approach is relatively robust even in configurations where the single cluster assumption is not entirely fulfilled, and (b) that the eigenvalue decomposition approach correctly identifies the simulated clusters even for low coupling strengths. Using the eigenvalue decomposition approach we studied spatiotemporal synchronization clusters in long-lasting multichannel EEG recordings from epilepsy patients and obtained results that fully confirm findings from well established neurophysiological examination techniques. Multivariate time series analysis methods such as synchronization cluster analysis that account for nonlinearities in the data are expected to provide complementary information which allows to gain deeper insights into the collective dynamics of spatially extended complex systems.


I. INTRODUCTION
Spatially extended complex dynamical systems may be thought of being composed of numerous constituents (dynamically formed subsystems) each having its own dynamics. Typically the relevant state variables of such systems can not be observed directly but only through some observation function that projects the high-dimensional state space onto an observation space of much lower dimension, resulting in a set of time series. Multivariate analyses of such time series might then help to gain deeper insights into the collective dynamics of spatially extended systems. Although a number of time series analysis methods have been developed over the past (see [1][2][3][4][5] for an overview) most techniques either allow to characterize single time series (univariate approaches) or to investigate relationships between two time series (bivariate approaches). However, applying bivariate techniques to pairs of time series -taken from a multichannel recording -does not necessarily allow to identify the relevant information in the full data set. The latter is of particular interest for scientific fields investigating spatially extended dynamical systems, such as meteorology, economics, social science or neurosciences, where a complex but relatively sparse connectivity between subsys-tems prevails. Understanding brain function -both during physiological and pathophysiological conditions (as e.g. in the case of epilepsy) -requires a characterization and quantification of the collective behavior of neural networks generating signals at different areas.
In principle, multivariate time series analysis techniques can be used to investigate mutual relationships between arbitrary numbers of time series. A large variety of methods [6] aim at revealing additional information by classifying time series into different groups. In addition to the classical principal component analysis (also known as Karhunen-Loeve transform) [7] independent component analysis [8] provides a decomposition of data into independent source signals, and if the assumption of independence holds, it can be regarded as a suitable method. If independence can not be assumed, mutual information based methods might be more appropriate [9,10]. The partial coherence [3] measures the fraction of coherence between two time series that is not shared with a third time series. Whereas the partial coherence is based on the assumption of linearity and thus does not capture nonlinear interactions, the recently proposed concept of partial phase synchronization [11] was designed to account for nonlinearities of the dynamics under investigation. In order to study causal relations among simultaneously acquired time series generated by linear stochastic systems the Granger causality [12] can be used by fitting autoregressive models. Besides a recently suggested nonlinear extension of Granger's ideas [13], we mention the directed transfer function that is defined for an arbitrary number of channels [14] and is based on a multivariate autoregressive model approach [15]. In Refs. [16,17] the partial directed coherence has been proposed for inference of linear Granger causality in the frequency domain based on vector autoregressive models of appropriate order. Both methods have been repeatedly applied to study interdependencies and causal relationships among neural signals (see e.g. [18][19][20][21] and references therein).
Over the last decade time series analysis techniques known from random matrix theory [22,23] have been repeatedly shown to allow an improved characterization of complex spatiotemporal correlation patterns. In these studies particularly the equal-time correlation matrix has been analyzed that was constructed from multivariate data sets obtained empirically in scientific fields ranging from economy [24,25] and meteorology [26] to the neurosciences [27][28][29][30][31]. More recently an approach for the detection of clusters in financial data based on properties of the eigenvalue spectrum of the equal-time correlation matrix was proposed in Ref. [32]. By filtering out the random part and collective market-wide effects the authors were able to detect groups of stocks by optimizing the matrix representation. Such an approach, however, requires a clear-cut definition of the random part, which is usually assumed to be associated with the small bulk eigenvalues. In Ref. [29] it was demonstrated though that the lower part of the eigenvalue spectrum may contain essential information.
Another way to study interactions in spatially extended systems is based on a statistical analysis of phase synchronization phenomena [4,33]. In most studies, however, the analysis of empirical multivariate data has been accomplished by a repeated application of bivariate synchronization measures and it remains to be established whether this approach allows to fully characterize a common integrating structure that may be present in the data. Addressing this issue, Allefeld and Kurths proposed a genuinely multivariate phase synchronization analysis method [34,35] and successfully applied their method to electroencephalographic (EEG) data recorded during a psychological experiment. The authors concluded their method to provide additional information on brain dynamics in a topographically, temporally, and frequency-specific way, as well as in other fields concerned with multivariate oscillatory processes. More recently, Allefeld and colleagues [36] introduced a new approach that addresses a limitation of their original method, namely the assumption of the existence of a single synchronization cluster in the data. By using methods from random matrix theory the resulting approach appears to be capable of identifying multiple synchronization clusters in the data which makes it highly attractive to characterize pathophysiological, spatiotemporal synchronization phenomena in multichannel EEG recordings from epilepsy patients.
Recent findings indicate that bivariate analysis techniques allow to characterize physiological and pathophys-iological phenomena in the human brain (see Ref. [37,38] for an overview), and it can be expected that multivariate phase synchronization analysis techniques provide complementary information. To address this issue we here study the synchronization cluster analysis methods proposed in Refs. [35,36] particularly with respect to field applications using model systems. In addition, we show that the method proposed in Ref. [36] allows to detect multiple synchronization clusters in long-lasting multichannel EEG recordings from epilepsy patients. This paper is organized as follows. Since both methods are based on the mean phase coherence we first recall its definition and interpretation as a bivariate measure for phase synchronization (Sec. II A). In Sec. II B we briefly introduce the multivariate synchronization cluster methods, namely the mean field approach (Sec. II B 1) and the eigenvalue decomposition approach (Sec. II B 2). Next, we present our simulation studies aiming at an exploration of the limitations and a comparison of both methods (Sec. III A). Finally, in Sec. III B we present findings that were obtained from a spatiotemporal synchronization cluster analysis of multichannel EEG data recorded from epilepsy patients using the eigenvalue decomposition approach.

A. Measuring phase synchronization
Phase synchronization was first described by Christiaan Huygens [39] in the 17th century and can be defined as the locking of the phases of two oscillating systems j and k: In a statistical way the degree of phase synchronization can be quantified by measuring the phase differences ∆ϕ jk n times and transforming them onto a unit cycle in the complex plane. The underlying circular distribution of the sample can be characterized by means of directional statistics [40] with the mean phase coherence R jk [41,42] where ϕ jm denotes the phase of system j in measurement m. By definition, R jk is confined to the interval [0,1] where R jk = 1 indicates fully synchronized systems. In field applications the sample size n is typically limited making R jk an estimate of the true population value ρ jk = e i(ϕj−ϕ k ) of the underlying distribution of phase differences (the angle brackets denote the average over all members of the population). When analyzing real valued time series s(t) different methods can be used to extract phase information. Methods based on the Fourier, the Hilbert, or the wavelet transform were shown to be equivalent under relatively general assumptions [43,44]. The main idea is to map the data onto the complex plane using a function z and to take the complex argument in order to obtain the phase ϕ = arg(z(t)). We here followed an approach based on the Hilbert transform using the analytic signal [45,46] where Hs(t) denotes the Hilbert transform of the signal s(t) and P denoting the Cauchy principal value of the integral. Application of the convolution theorem turns the last equation into where F denotes the Fourier transformation and F −1 the inverse transformation, respectively. Thus, the imaginary part of the analytic signal is obtained by shifting each frequency component of the original signal by π/2 separately. It is important to keep in mind that phases may not have a physical meaning for arbitrary signals. If, however, the dynamics exhibit oscillations with a single main rhythm, then the phases are typically well defined [47,48]. In our applications (see Sec. III) we computed R jk from the discrete time series s j (t) and s k (t) as follows. First, the data were normalized to zero mean, which corresponds to setting the DC Fourier coefficient to zero. Second, in order to avoid edge effects the first and last 10% of a data window of size n ′ were tapered using a cosine half wave (Hanning window) before performing the Fourier transform. Third, since the computation of the Hilbert transform requires integration over infinite time, which cannot be performed for a window of finite length, 10% of the calculated phase time series ϕ j (t) and ϕ k (t) were discarded on each side of the window, reducing the size of the considered phase time series to n = 0.8 · n ′ . For a data window of size n and m ∈ {0, . . . , n − 1} denoting the data point within the window the mean phase coherence R jk between the two time series with sampling interval ∆t was obtained by identifying ϕ jm = ϕ j (m∆t) and applying Eq. (2).

B. Synchronization cluster analysis
For N oscillating systems the pairwise computation of the mean phase coherence R jk leads to a matrix R which is symmetric due to definition (2). Subsets of oscillating systems can be interpreted as synchronization clusters if these systems exhibit higher mean phase coherence values between each other than with systems not included in the same subset. Allefeld and colleagues [35,36] recently proposed two different multivariate approaches that aim at identifying synchronization clusters in which the oscillating systems participate with different strength. Both approaches are based on the bivariate mean phase coherence and shall briefly be recalled here.

Mean field approach
In Ref. [35] a mean field approach has been presented which assumes the existence of a single synchronization cluster C in the data. Generating a common rhythm all oscillating systems constitute the cluster but contribute to its emergence to a different extent. A mean field with phase Φ can then characterize the dynamics of the collective behavior. Using Eq. (2) the degree of participation (participation strength) of each oscillating system j to the cluster C can be quantified by A straightforward derivation of the phase Φ has been demonstrated in Ref. [35] for a simple model system. This method, however, requires exact knowledge about the underlying equations of motion and thus can not be applied to unknown systems in general. Nevertheless, the participation strength (6) can be estimated directly as follows.
If the assumption of the existence of a single cluster holds, a mean field can be introduced such that the dynamics of the phase differences are decoupled. If, in addition, the noise affecting the oscillating systems is statistically independent for each system the phase differences ∆ϕ j = ϕ j − Φ become independent random variables. Hence the population values ρ jk of the mean phase coherences R jk turn into Taking into account that the mean phase coherence is asymptotically normally distributed R jk ∼ N (ρ jk , σ jk ) [49] and is an empirical estimate of ρ jk = ρ jC ρ kC (see Sect. II A), a maximum likelihood estimation of the ρ jC leads to the minimization of the sum of squared weighted errors that defines the cost function Assuming that the circular distribution of the phase differences can be sufficiently approximated by a wrapped normal distribution, the standard deviation σ jk of the sampling distribution of R jk can be expressed in terms of ρ jk (cf. [40,49]) In the following we refer to ρ jC as defined above as R jC , thereby emphasizing an interpretation as a to-cluster synchronization strength analogous to Eq. (6). In our applications participation strengths were computed by minimizing Γ using an iterative algorithm proposed in Ref. [50].

Eigenvalue decomposition approach
In Ref. [36] another approach has recently been proposed that is based on the eigenvalue decomposition of the matrix R and appears to allow identification of multiple clusters. The procedure makes use of findings by Müller and colleagues [29] who demonstrated that information about the correlation structure of multivariate data sets is imprinted into the dynamics of the eigenvalues and into the structure of the corresponding eigenvectors by nonrandom level repulsion.
The eigenvalues λ c and eigenvectors ν c of the symmetric and real-valued matrix R (see Eq. 2) can be obtained by solving the eigenvalue equation which, in general, has N different solutions. In the following it is assumed, that the eigenvectors are normalized (| ν c | = 1). Being transformed by an orthogonal transformation into the basis of its eigenvectors R becomes the diagonal matrix D of its eigenvalues. The invariance of the trace under this transformation leads to the equation N = tr(R) = tr(D) = c λ c . In the case of systems j, k showing no phase synchronization (R jk = 0 for j = k) the equation is trivially fulfilled by λ c = 1 for all c, whereas the occurrence of entries R jk > 0 for j = k induces a level repulsion, a combined increase and decrease of eigenvalues in such a way that N = c λ c still holds. A multivariate analysis is realized as follows (cf. Ref. [36]): Each eigenvalue λ c > 1 is associated with a synchronization cluster and quantifies its strength within the data set. The internal structure of cluster c is described by the corresponding eigenvector ν c . Being normalized ( j ν 2 jc = 1) its components quantify the relative involvement of each system j to cluster c by ν 2 jc . Combining the eigenvalue λ c and the index ν 2 jc the "absolute" involvement of a system j in a cluster c can be described by the participation index Consequently, system j is considered as belonging to cluster c for which its participation index becomes maximal.
In case of non-vanishing entries of R between systems belonging to different clusters (inter-cluster synchronization) it can be observed that these clusters are not characterized by separate but by a superposition of eigenvectors. In order to adjust the interpretation of the participation indices as the degree of involvement of a system within one cluster (see. Eq. 11) it was proposed in Ref.
[36] to compute the p j,c in a first step and to assign the systems to the clusters as mentioned above. In a second step the matrix entries representing inter-cluster synchronization are set to zero and the participation indices are computed on the trimmed matrix again. In our applications we followed exactly this scheme.
Summarizing this section we conclude that both methods, the mean field and the eigenvalue decomposition approach, seem to provide similar results in single cluster configurations for which an almost functional dependency of R 2 jC = p j,1 (λ 1 denoting the largest eigenvalue) can be observed [36]. Nevertheless, the question whether the approaches provide meaningful results in case of multi-cluster configurations -particularly with regard to field applications -remains unaddressed and shall be investigated in the following section.

A. Simulations
We here studied a three-cluster configuration which consisted of a lattice of 32 coupled identical Lorenz systems (cf. Fig. 1). Each system j is defined by the differential equationṡ where ǫ j denotes the coupling strength controlling the influence of the diffusive unidirectional coupling on system j by one of three chosen driving systems D 1,2,3 . In order to control this configuration by a single parameter we set ǫ j = ǫ for systems within each cluster while ǫ j = 0 for the uncoupled systems. Taking randomly chosen points in the state space near the Lorenz attractor as initial conditions the differential equations were iterated using a fourth order Runge-Kutta algorithm [51] with a step size of 0.01. In order to eliminate transients, the first 10 4 iterations were discarded. For increasing coupling  strengths ǫ = 0.0 to 1.4 (step size 0.05) we generated scalar time series of the x-components (n ′ = 5 · 10 5 data points) and computed the participation strength R jC and the participation indices p j,c .
For the mean field approach the configuration represents a violation of the single-cluster assumption. When increasing the coupling strength ǫ the mean phase coherence values between the coupled systems increased and consequently the three synchronization clusters emerged. This is reflected by the dependency of the cost function Γ on ǫ as shown in Fig. 2a. Since the single-cluster assumption leading to Eq. (7) was violated, Γ increased rapidly for ǫ > 0.4. For comparison, we repeated the analysis for a single-cluster configuration that was generated by setting the coupling terms between the systems of clusters driven by D 1 and D 2 to zero. Here Γ ≤ 0.6 for all ǫ values (cf. Fig. 2a). Using Γ as an indicator that reflects the violation of the single-cluster assumption (cf. Ref. [36]), which increases with the internal degree of synchronization between systems within a cluster, we generated plots (cf. Fig. 2b) showing the spatial distribution of the participation strength R jC for the three-cluster configuration chosen here. For ǫ ≤ 0.4 (Γ ≤ 0.6) no clear-cut cluster structure could be identified in the spatial distribution of participation strengths. For ǫ ≥ 0.5 (Γ ≫ 0.6) either a single or all three synchronization clusters emerged, however, with a varying degree of visibility. Interestingly, for various coupling strengths (Fig. 2b, e.g. ǫ ∈ {0.6, 0.9, 1.1}) the participation strengths of systems involved in one cluster exhibited higher values than the remaining participation strengths. This behavior is shown in more detail in Fig.  3a for exemplary systems involved in one cluster. In order to further elucidate this phenomenon we calculated the sum S d = j>k;j,k∈C d R jk , where C d denotes the index set of systems belonging to cluster d. S d quantifies the strength of cluster d and increased, on average, with an increasing coupling strength ǫ (cf. Fig. 3b) while fluctuations can be attributed to the mean phase coherence values being asymptotically normally distributed. A synchronization cluster d became visible in the spatial distribution of participation strengths R jC when its strength S d dominated the configuration (cf. Fig. 3). This effect was caused by the cost function Γ exhibiting several local minima whose values changed according to the cluster strengths S d . For ǫ = 0.8 all S d attained similar values. Here the minimization algorithm terminated because of having found a single minimum which would lead to the (mis)interpretation of all systems belonging to one single cluster.
In contrast to the mean field method, the eigenvalue decomposition approach provides not just a single set of scalars but multiple sets of scalars representing different clusters. The eigenvalues were labeled here according to λ 1 ≥ . . . ≥ λ N thereby enabling an easy identification of the dominating cluster structures by the corresponding index c. As shown in Fig. 4 for low coupling strengths (ǫ ≤ 0.3) the method was not able to detect the three different clusters but indicated a multitude of small clusters. This can be attributed to random fluctuations of the mean phase coherences. Increasing slightly the coupling strength to ǫ = 0.4 two of the three simulated clusters were already detected in the two largest eigenvalues λ 1 and λ 2 whereas the corresponding participation strengths did not reveal any cluster structure (cf. Fig. 2b). For ǫ ≥ 0.5 all three clusters were clearly visible and identified by the three largest eigenvalues. Their participation indices increased with increasing coupling strength ǫ (see e.g. ǫ = 1.1 to 1.4) thereby reflecting the internal degree of synchronization within the clusters. The cluster indices c changed for different coupling strengths ǫ due to the fact that the cluster strengths fluctuated (cf. S d in Fig. 3b) which is reflected by the eigenvalues (see e.g. ǫ ∈ {1.1, 1.3} in Fig. 3b and Fig. 4).
In contrast to the mean field approach the eigenvalue decomposition approach was capable of distinguishing between different clusters (see [36]). Limitations of the method, however, that are related to parameters influencing the detectability of different clusters have not been studied so far. To address this issue we designed a configuration consisting of N = 32 systems that formed two clusters C 1 = {1, . . . , r} and C 2 = {r + 1, . . . , N } where the parameter r controls the relative size of an individual cluster. We did not study dynamical systems here but instead assigned values to the entries of matrix R that were drawn from normal distributions given that the mean phase coherences are asymptotically normally distributed (see section II B 1). Mean phase coherences between systems of different clusters were drawn from the inter-cluster distribution N (ρ (int) , σ (int) ), while those representing the entries between systems of the same cluster C 1 and C 2 were drawn from N (ρ (1) , σ (1) ) and N (ρ (2) , σ (2) ), respectively. For a given population value ρ the standard deviation σ was determined by equation (9) using a sample size of n = 200.
Recalling that the detection of different clusters is based on the participation indices computed in the first processing step of the method (see Sec. II B 2), which assigns a system to a cluster for which its participation index becomes maximal, we here quantify an erroneous assignment by considering the differences of the participation indices obtained from the first processing step χ j12 = p j,1 − p j,2 of system j and clusters C 1 and C 2 . For the configuration considered here we define the weighted amount of erroneous assignments as ∆ = ∆ 1 + ∆ 2 : N λ = 2 0 : else (13) where N λ denotes the number of eigenvalues λ > 1 and ∆ 1 and ∆ 2 are determined by where Θ is the Heaviside step function. If system j is erroneously assigned to the cluster C 1 but belongs to a cluster C 2 by construction, Θ(χ j12 ) will be larger than zero. Thus ∆ 2 (∆ 1 ) sums up incorrect assignments of systems j to cluster C 1 (C 2 ) weighted by the differences of the corresponding participation indices χ. If all systems are assigned to the correct cluster, then ∆ = 0 by definition. Since the eigenvalues quantify the cluster strengths, a variation of the relative size of an individual cluster (by varying r) makes the largest eigenvalue traverse a minimum. This enables us to label the eigenvalues λ c > 1, c ∈ {1, 2} according to their clusters C 1 or C 2 . We studied a transition from a configuration that consisted of two clusters C 1 and C 2 (ρ (1) = ρ (2) = 0.8) to a single-cluster configuration C 0 = {1, . . . , N } by successively increasing ρ (int) from 0.0 to 0.8 in steps of 0.01. Fig. 5 shows the weighted amount of erroneous assignments ∆ depending on the relative size of an individual cluster r and on the inter-cluster synchronization level ρ (int) . The eigenvalue decomposition method successfully identified (∆ = 0) the two clusters for low values of the inter-cluster synchronization level (ρ (int) < 0.17) except for clusters of equal size (r = N/2). The weighted amount of erroneous assignments rapidly increased, however, for increasing ρ (int) , thereby failing to detect the  cluster structures for an increasing set of r values. The number of eigenvalues N λ being larger than 1 appeared to be a sensitive measure for the number of clusters present in the data. We observed N λ → 1 (i.e., ∆ = 0 by definition, see Eq. 13) only for relatively large values of the inter-cluster synchronization level (ρ (int) > 0.7) where the transition to a single-cluster configuration was completed.
Next we studied a transition from a configuration of two clusters C 1 and C 2 (ρ (1) = ρ (2) = 0.8, ρ (int) = 0.2) to a configuration with only one cluster (here C 1 ) by decreasing ρ (2) in steps of 0.01 down to 0.2 (see Fig. 6a). The weighted amount of erroneous assignments ∆ vanished for ρ (2) approaching the mean synchronization level of the inter-cluster distribution (ρ (int) = 0.2). In this range the number of eigenvalues N λ increased (Fig. 6b). These additional structures were formed by those systems, which were previously involved in cluster C 2 for higher values of ρ (2) . This fragmentation of a cluster c can also be observed when increasing the corresponding standard deviation σ (c) independently from the population value ρ (c) , thereby taking into account various uncertainties (e.g. definition of observables, measurement precision, finite sample size, phase extraction). Since the additional structures were caused by random fluctuations of the mean phase coherences within the cluster, we refer to them as pseudo-clusters in the following. In our simulations we observed pseudo-clusters of different sizes, whose corresponding eigenvalues were only slightly larger than unity compared to the eigenvalues of the true clusters.
When varying ρ (2) the set of r values for which ∆ > 0 remained almost stable but was shifted to lower r values (Fig. 6a). Therefore, together with the relative size of an individual cluster the level of synchronization between systems within a cluster must be regarded as crucial for the eigenvalue decomposition approach. The results shown in Figs. 5 and 6a can then be sufficiently explained with the help of the cluster strength S d which considers both, the cluster size and the level of synchronization between systems within a cluster. The method failed to correctly assign all systems to their clusters in case of similar cluster strengths |S 1 − S 2 | < κ, where κ increases with an increasing level of inter-cluster synchronization. The latter leads to superpositions in the eigenvector components and, as a result, makes it impossible for the method to correctly detect the clusters by the maximum participation index criterion in certain configurations. In order to further elucidate this phenomenon we show in Fig. 7 the components ν jc of the eigenvectors corresponding to the two largest eigenvalues for the case r = 16 and ρ (2) = 0.8. The squared components ν 2 jc of both vectors attained similar values, causing the method to classify all systems as belonging to a single cluster. Nevertheless, the two clusters C 1 and C 2 were clearly visible in the linear combinations ν 1 ± ν 2 .
Summarizing this section we conclude that the eigenvalue decomposition approach successfully identified the three synchronization clusters in our simulations using coupled Lorenz systems and, moreover, appears to be sensitive to detect clusters even for low coupling strengths. Although the mean field approach is by definition not able to distinguish between different clusters, the method appears to be relatively robust even in situations where the single-cluster assumption is not (entirely) fulfilled, leading to higher participation strengths for sys- tems of the dominating cluster in our simulation. When comparing both methods the eigenvalue decomposition approach can be regarded as superior when analyzing data exhibiting different synchronization clusters, which can be expected e.g. for multichannel EEG data.

B. EEG data
In this section we present exemplary findings obtained from applying the eigenvalue decomposition method to multichannel EEG time series that were recorded from three epilepsy patients (denoted as P1, P2, and P3) suffering from pharmacoresistant focal epilepsies of neocortical origin. For these patients complete seizure control can be obtained by resecting the part of the brain responsible for seizure generation (epileptic focus). This requires an exact localization of the epileptic focus and its delineation from functionally relevant brain structures during the presurgical workup. When no concordant information can be achieved from noninvasive diagnostic techniques, the EEG is recorded from implanted electrodes over a longer period, typically 2-3 weeks. The analyses reported here were made after surgery had taken place, and after it had become clear from its success whether the localization of the epileptic focus had been correctly predicted. All patients had signed informed consent that their clinical data might be used and published for research purposes, and the study was approved by the local medical ethics committee.
Previous studies have shown that even during seizurefree intervals the seizure generating area of the brain exhibited higher interdependences [52] and an higher degree of synchronization [42] than other brain areas. Together with results obtained from applying univariate time series analysis techniques (see e.g. [53] and references therein) these findings allow an improved understanding of intermittent dysfunctioning of the brain between seizures and provide potentially useful diagnostic information. We here addressed the question whether complementary information can also be obtained from a multivariate approach. Specifically, we investigated whether the eigenvalue decomposition approach (a) allows to localize the epileptic focus analyzing EEG recordings from the seizure-free interval only, and (b) is capable of detecting short time changes of synchronization patterns associated with physiological processes in the human brain. To this end we analyzed continuous EEG recordings that lasted 39 h for patient P1 and 26 h for patient P2 covering different physiological and pathophysiological states of the patients. In addition, we analyzed an EEG recording lasting 30 minutes during which patient P3 was simultaneously presented item pairs, either word pairs or unpronounceable letter string pairs of 5-11 letters length, and was instructed to press a button if he/she recognized word pairs with identical or highly similar meanings (synonyms) (task T1) or identical letter strings (task T2). No button press was demanded in case of non-synonym word pairs or letter strings in which one consonant differed. Tasks alternated every 3 min. The whole sequence was adopted from Ref. [54] and comprised three blocks of word and letter string pairs each. A baseline recording of 10 minutes was performed with eyes open before and after the experiment. Recent findings [55,56] indicate that tasks such as T2 are associated with a greater demand for decision-making processes due to the involvement of different phenomena like reading, phonological retrieval, or orthographic analysis during such conditions. We thus hypothesized to observe a higher degree of collectivity among neuronal assemblies particularly in brain areas associated with language processing. We did not expect to observe synchronization phenomena associated with the epileptogenic process in these brain structures from patient P3 since the epileptic focus was localized in a more distant brain structure.
The EEG was measured from grid electrodes (rectangular flexible grids of N = 8 × 4 contacts) placed onto the temporal lateral neocortex (see Fig. 8). EEG data were sampled at 200 Hz using a 16-bit analog-to-digital converter and filtered within a frequency band of 0.3-70 Hz. Using a moving-window technique EEG signals were divided into segments of n ′ = 4096 sampling points each, and segments overlapped by 20%. The length of the re- sulting segments corresponded to 16.38 s at the given sampling rate and can be regarded as a compromise between the required statistical accuracy for the calculation of the mean phase coherence and the approximate stationarity within a segments length (see Ref. [57] and references therein). After the calculation of the mean phase coherence (cf. Sect. II A) for all channel combinations from each segment, the eigenvalue decomposition method was applied (cf. Sect. II B 2).
In order to estimate the number of clusters being detected on average in the data we sorted eigenvalues and corresponding eigenvectors of each segment in a descending order (λ 1 ≥ . . . ≥ λ N ) and computed the mean valuesλ c over all segments for each patient. The five to six largest averaged eigenvaluesλ c were larger than 1 causing the method to detect on average an according number of clusters per segment. Based on our findings presented in Sect. III A we expected the lower part of the eigenvalue spectrum being larger than 1 to yield a considerable amount of pseudo-clusters. We therefore restricted further analyses to the three largest eigenvalues.
Since a cluster is represented by a set of participating channels, each cluster can be written as a bit string of size N where the bits represent the channels being involved (1) or not involved (0) in the cluster. This allowed us to discard the information about the relative involvement of each system within one cluster as reflected by the participation indices and to handle the clusters in a convenient way, namely by considering their simplified representation. Bit strings connected to eigenvalues fulfilling the threshold criterion (λ > 1) varied largely over time as shown exemplarily in Fig. 9. Apart from frequently appearing bit strings in which all channels were involved, structures in time could be observed, where successive bit strings differed only in few bits. Given that the mean phase coherence fluctuates (due to measurement noise, a limited number of data points, or even due to short-term physiological or pathophysiological phenomena within a segment), the same cluster structure cannot be expected to show up in exactly the same bit string representation but may vary in some bits. In order to minimize this ef- fect we sorted bit strings that were visible for more than 6 minutes in the data into groups which differed only in up to three bits. Moreover, we discarded groups of bit strings which represented clusters of size 0, 1, or 32 channels. In Fig. 10 we show groups that represented the largest number of bit strings in the EEG data. Interestingly, for both patients the spatial distribution of electrodes involved in the cluster labeled "a" corresponded to the spatial extent of the epileptic focus as determined by the presurgical workup. This cluster could be observed throughout the datasets and in total for at least 112 minutes for patient P1 and 125 minutes for patient P2. Indeed, the patients were operated on exactly this region and are now free of seizures. Thus, we conclude that the method seems to be quite sensitive for detecting synchronization clusters in EEG time series recorded from epilepsy patients even during the seizure-free interval.
Whereas the aforementioned results were obtained by analyzing the bit string groups which corresponds to a temporal average (over all segments), we now investigated whether the eigenvalue decomposition approach allows to detect short-term changes of cluster structures that can be related to physiological synchronization phenomena (the language processing paradigm mentioned above). The participation indices were computed from the EEG recording from patient P3 and translated into bit strings using the same criteria as mentioned above. Bit strings corresponding to the three largest eigenvalues are shown in Fig. 11. Interestingly, a modulation of the bit strings depending on the tasks can be observed. The occurrence of bit strings representing a cluster of the electrode contacts A7, B7-8, and C7-8 is noticeable particularly during task T2 (letter matching task). For this patient the brain area covered by these electrode contacts was associated with language processing (Wernicke's area). The observed higher level of synchronization probably reflects the higher degree of collectivity among neuronal assemblies that are involved in the cognitive operations comprising this task. This suggests that the eigenvalue decomposition method is capable of detecting short time changes of synchronization patterns associated with physiological processes in the human brain.  11: (Color online) Spatio-temporal evolution of clusters defined by the largest, second-, and third-largest eigenvalues calculated from EEG data from patient P3 recorded during a language processing paradigm involving two different tasks (separated by black lines). The patient was simultaneously presented item pairs, either word pairs or unpronounceable letter string pairs, and was instructed to press a button upon recognizing synonyms (T1) or identical letter strings (T2). A baseline recording (B) was performed before and after the experiment. Each pixel in the figure represents the color-coded cluster membership of a given electrode contact (A1-A8, B1-B8, C1-C8, D1-D8 shown on the ordinate; cf. Fig. 8) and for a given EEG segment of duration 16.38 s (abscissa). Contacts A7, B7-8, and C7-8 covered a brain area associated with language processing (Wernicke's area) and were most noticeable during task T2.

IV. CONCLUSION
In this paper we have studied two multivariate phase synchronization analysis methods, namely a mean field approach [35] and an eigenvalue decomposition approach [36]. While the mean field approach assumes the existence of a single synchronization cluster in the data, the eigenvalue decomposition approach appears to be capable of identifying multiple clusters. Based on the results of numerical simulations of multiple synchronization clusters within a lattice of 32 coupled identical Lorenz systems, we demonstrated that the mean field approach appears to be relatively robust even in situations where the single-cluster assumption is not entirely fulfilled. The eigenvalue decomposition approach successfully identified multiple synchronization clusters in our simulations and appears to be sensitive to detect clusters even for weak couplings. However, in case of non-vanishing intercluster synchronization the method failed to correctly assign the systems to their clusters in certain configurations. The influence of measurement noise, which was not discussed in depth here, can be expected to lead to the fragmentation of clusters present in the data. Nevertheless, the eigenvalue decomposition approach can be regarded as superior when analyzing data exhibiting different synchronization clusters.
When being applied to field data a number of influenc-ing factors limit the significance of the eigenvalue decomposition approach and need further investigations. These factors include the finite length of available data, random spatio-temporal correlations, or even non-random correlations being induced by the data acquisition system (e.g. filtering or, as in the case of EEG recordings, the choice of a suitable reference). Nevertheless, our preliminary applications of the eigenvalue decomposition approach to multichannel EEG recordings from epilepsy patients indicate that the method allows to gain deeper insights into the collective dynamics of neuronal networks, both under physiological and pathophysiological conditions. Despite the limited EEG database used in this study, the achieved results can be regarded as promising. Further evaluations on a larger EEG database are currently underway.