## Abstract

The purpose of the present paper is to develop a method, based on equal-time correlation, correlation matrix analysis and surrogate resampling, that is able to quantify and describe properties of synchronization of population neuronal activity recorded simultaneously from multiple sites. Initially, Lorenz-type oscillators were used to model multiple time series with different patterns of synchronization. Eigenvalue and eigenvector decomposition was then applied to identify “clusters” of locally synchronized activity and to calculate a “global synchronization index.” This method was then applied to multichannel data recorded from an in vitro model of epileptic seizures. The results demonstrate that this novel method can be successfully used to analyze synchronization between multiple neuronal population series.

## INTRODUCTION

Neuronal population synchronization is an important feature of normal and abnormal brain function (Buzsáki 2004; Buzsáki and Draguhn 2004; Engel et al. 2001; Rosso et al. 2001; Traub et al. 1998, 1999). Epileptic seizures, in particular, are associated with excessively synchronized rhythmic discharges in large populations of neurons (Jefferys and Fox 2007) and it has been suggested that, immediately preceding a seizure, there is a fusion of activity initially generated in small independent clusters of neurons (Bikson et al. 2003; Fox et al. 2007). Further understanding of the mechanisms underlying synchronization between neuronal populations would be facilitated by the development of computing techniques that can estimate the degree of synchronization within and between neuronal populations of varying size (Le Van Quyen and Bragin 2007).

Existing methods to measure synchronization of neuronal population activity include cross-correlation, mutual information, spectrum-based coherence, nonlinear interdependence, and phase synchronization (Brown and Kocarev 2000; Carmeli et al. 2005; Le Van Quyen et al. 2005; Mormann et al. 2000; Varela et al. 2007). These methods, however, cannot be applied to more than two series. Recently, methods based on random matrix theory (RMT) have been developed to analyze multiple time series (Allefeld et al. 2007; Müller and Baier 2005; Müller et al. 2006). The first application was to investigate the statistics of energy levels of complex quantum systems (Brody et al. 1981; Guhr et al. 1998; Mehta 1991), but it has subsequently been used to analyze financial data (Plerou et al. 1999; Wilcox and Gebbie 2007), magnetoencephalographic recordings (Kwapien et al. 1998, 2000; Tass et al. 2003), electroencephalographic recordings (Allefeld et al. 2007; Müller and Baier 2005; Müller et al. 2006; Seba 2003), climate data (Mayya and Amritkar 2006; Santhanam and Patra 2001), array signals (Everson and Roberts 2000), and networks (Goh et al. 2001; Luo et al. 2006).

In this report, we use an equal-time correlation method to estimate correlation between all pairs of signals to produce a correlation matrix from *1*) multiple Lorenz-type oscillators and *2*) neuronal population activity simultaneously recorded with multiple electrodes. Then, we use a method based on RMT to analyze the matrices, which we term *correlation matrix analysis* (CMA): we demonstrate that this novel method can successfully extract detailed information about cluster formation, number and composition of clusters, global synchronization, and how synchronization patterns change with time. The present analysis is based on field potentials that are produced by synchronous activity in populations of neurons close to each recording electrode. Other methods have been successfully used to investigate synchrony across large networks of spiking neurons (e.g., Ginzburg and Sompolinsky 1994; Golomb 2007; Hansel and Sompolinsky 1996); in principle, multiple spike trains generated by single cells could also be analyzed using CMA, but the choice of appropriate correlation methods is not explored here.

## METHODS

### Equal-time correlation and CMA methods

Selection of the most appropriate method, for construction of the correlation matrix, depends on signal characteristics; when the time delay between two neural population activity series is large, phase synchronization is appropriate; when the time delay between two series is very small (e.g., neural activity series recorded in brain tissue slices in vitro), equal-time correlation can be applied. The equal-time correlation is a simple method for measuring synchronization. The multivariate neuronal activity series is *Z _{i}*(

*t*), where

_{k}*i*= 1, …,

*M*over a time window

*t*(

_{k}*k*= 1, …,

*T*);

*M*is the channel number; and

*T*is the number of data points in the time window (

*T*should be considerably larger than

*M*). To provide the same scale for all of neuronal population activity, the data set is normalized by (1) where 〈

*Z*〉 and σ

_{i}_{i}are the mean and the SD of

*Z*(

_{i}*t*), respectively. An equal-time correlation matrix C can then be constructed as (2) Due to the normalization all entries of matrix

_{k}**C**range from −1 to 1: when

*c*= 1, there is a perfect correlation between the

_{ij}*i*and

*j*series; when

*c*= −1, there is a perfect anticorrelation; and when

_{ij}*c*= 0, there is no correlation. Thus all diagonal elements are equal to 1 because each signal is perfectly autocorrelated. The eigenvalue decomposition of

_{ij}**C**is given by (3) where λ

_{i}is the eigenvalue and λ

_{1}≤ λ

_{2}≤… ≤ λ

_{M};

**v**

_{i}is the eigenvector corresponding to λ

_{i}. The eigenvalues have the following properties:

*1*) All eigenvalues are real numbers and the sum of the eigenvalues equals the number of series

*M*because

**C**is a real symmetric matrix (e.g., for a nine-channel recording, the sum of eigenvalues = 9).

*2*) If the time series are fully uncorrelated (nonsynchronized),

**C**will approximate to an identical matrix and all of the eigenvalues will distribute around the value 1 (indicating the presence of random correlations between the time series).

*3*) Once all of time series are completely correlated, the elements of

**C**are equal to 1 and the maximum eigenvalue is equal to the number of time series

*M*, other eigenvalues falling to zero. Thus eigenvalues can provide information about the synchronization between individual elements of the matrix and the highest eigenvalue can give information about global synchronization.

To obtain a normalized value of global synchrony (which would range from 0 to 1) and be independent of the number of channels (*M*) we derived a global synchronization index as follows. First, we randomize all channels of the time series, thereby destroying all the equal-time correlations that exist; then compute a surrogate correlation matrix **R**; the eigenvalues of this matrix **R** can be obtained from *λ _{k}* (

*k*= 1, …,

*M*). Repeating this randomization by

*N*, the mean and SD of eigenvalues are

*λ̄*and

_{k}*SD*. The randomization can be completed by a resampling technique; different correlation methods may need different resampling methods (Kantz and Schreiber 2003). In the present work, we randomize the phase relationship using the amplitude-adjusted Fourier transform (AAFT) (Schreiber and Schmitz 1996). The original series is denoted by the vector {

_{k}*y*(

*n*)}. First, the Fourier spectrum of this vector is calculated and a phase series is obtained. To preserve the power spectrum in the surrogate series, the squared coefficients of the surrogate series are replaced by those of the original time series, the phases in the Fourier spectrum are randomized, and, finally, an inverse Fourier transform is applied to construct a new series. The AAFT algorithm is subsequently explained mathematically. The linear properties of the original series

*y*(

*n*) can be specified by the squared amplitude of its discrete Fourier transform (4) The surrogate time series

*y*′(

*n*) can be created by multiplying the Fourier transform of the original series by random phases and then inverse Fourier transform to the time domain (5) where 0 ≤ θ

_{k}< 2π represents independent uniform random numbers. Thus the new series will have the same power spectrum as the original series, but the phase is different.

Once we obtain the *λ̄ _{k}*, the normalized synchronization can be computed by the following equation (6) where

*k*= 1, …,

*M*. The synchronization index,

*Syn*_

*Index*, ranges from 0 to 1, where 0 denotes no synchrony and 1 denotes perfect synchrony among the time series. The maximal value of the synchronization index is called the global synchronization index.

_{k}*K*is a constant that determines the threshold. To allow for the nine comparisons used in this analysis, a Bonferroni correction (SISA 2007) gives a

*K*value of 3 for an overall significance of 0.01.

In addition to the synchronization index, the eigenvalues can give information about the number of clusters of synchronization. Using the resampling method, the surrogate eigenvalues and SD of the random time series are computed. If the original eigenvalue (λ_{k}) is greater than the sum of the averaged surrogate eigenvalues (*λ̄ _{k}*) and SD (

*SD*) ×

_{k}*K*(i.e.,

*λ̄*+

_{k}*K*×

*SD*), then it suggests the presence of clusters of synchronization between at least two elements of the synchronization matrix. The equation to compute the number of clusters is expressed as (7) where

_{k}*sgn*is a sign function [i.e., if λ

_{k}> (

*λ̄*+

_{k}*K*×

*SD*) is true, sgn is 1; otherwise, it is 0].

_{k}The eigenvalue can indicate the strength of a cluster. Corresponding to the eigenvalue, the eigenvector can describe its internal structure and *∑ _{i} v_{ik}^{2}*. In fact, the index

*v*is the weight value of each channel

_{ik}^{2}*i*involved to cluster

*k*. Thus qualification of the structure of synchronized clusters should be described by combining the eigenvalues and eigenvectors, in a variable called the

*participation index*(PI; Allefeld et al. 2007), which is defined as (8) where

*v*is the

_{ik}*i*th element of eigenvector

*v*and λ

_{k}_{k}is a corresponding eigenvalue. In other words, PI

_{ik}describes the contribution of the

*i*th series to the cluster

*k*. In terms of the distribution of the PI

_{ik}, the element of the cluster can be indicated, as can be seen in Fig. 1

*C*; in turn it is known which channels contribute to the cluster.

### Simulation analysis

To test the properties of the algorithm, initially a chaos system of coupled Lorenz-type oscillators (Veeramani et al. 2004) was simulated to generate a multivariate time series. The equation of the oscillator *i* in the Lorenz system is (9) where *i* is the index for the oscillator (*i* = 1, …, *M*; *M* = 10, which is the number of oscillators in the dynamical system). σ, γ, β are the parameters of the oscillators; (σ, γ, β) = (175, 10, 8/3) in this study. ε_{i}_{,j} (*i*, *j* = 1, …, 10) denotes the coupling strengths between oscillator *j* and oscillator *i*; the coupling coefficient for correlation or noncorrelation is either ε_{i}_{,j} = 1 or ε_{i}_{,j} = 0. The length of all ten series is 100 points. By changing the coupling coefficients ε_{i}_{,j}, four cases are considered: *1*) all the time series are random, i.e., ε_{i}_{,j} = 0; thus it models that no relationship between time series exists; *2*) time series 1–4 are correlated, i.e., ε_{i}_{,j} = 1 (*i*, *j* ≤ 4); and the others are random, i.e., ε_{i}_{,j} = 0 (*i*, *j* > 4); *3*) the time series 1–4 form one synchronization cluster and the time series 5–7 form a second cluster; and *4*) all the time series are correlated, i.e., ε_{i}_{,j} = 1 (*i*, *j* = 1, …, 10). The reason for selecting the Lorenz system is that it can generate oscillatory time series that can be coupled selectively. The system thus generates different correlation structures that are useful for testing our analytical methods.

Figure 1 shows the correlation structure of four case studies, including the equal-time correlation matrices, the eigenvalues, and the participation index. As can be seen from Fig. 1*A*, these four cases have different correlation patterns. Figure 1*B* plots the eigenvalues for the four cases; the number of eigenvalues for each case is equal to the number of channels and they are listed in order from λ_{1} to λ_{10}, where λ_{1} is the smallest and λ_{10} is the largest. In the random case (case 1), the eigenvalues increase gradually (∼1) and in order; none of the eigenvalues is significantly larger than the others (none of them is above the threshold of averaged surrogate eigenvalues) and the deviation of eigenvalues is small. By contrast, in case 2 there is one large eigenvalue (the maximal eigenvalue is λ_{10}); this is derived from the correlation of channels 1 to 4 (note that there is a proportionate reduction in eigenvalues λ_{1} to λ_{3}). Case 3 reveals two significant (above threshold) eigenvalues (λ_{10} and λ_{9}), which are derived from two independent clusters (i.e., the number of eigenvalues that cross the threshold value indicates the number of independent clusters of synchronized activity). Case 4 shows a complete correlation pattern; here the largest eigenvalue (λ_{10}) approaches 10 and the other eigenvalues are very small.

Eigenvalues can indicate the number of clusters and corresponding synchronization strength. The next point of interest is to determine the structure of individual clusters; for this, we use the participation index (PI). In Fig. 1*C*, case 1, there are ten curves; each eigenvalue corresponds to one curve. PI(10) shows the participation index of the maximal eigenvalue; the *x*-axis is the number of channel. The finding is that the PIs are randomly distributed: this is because the multivariate series are independent. In case 2, eigenvalues demonstrate the presence of one cluster; the participation indices PI(10) show that this cluster results from synchronization between series 1 to 4 [PI(10) corresponds to the largest eigenvalue λ_{10}]. This is in agreement with the settings that were used to simulate this multivariate time series. In case 3, we can find two clusters, the participation indices PI(10) indicate one cluster that is composed of series 1 to 4 (corresponding to eigenvalue λ_{10}) and indices PI(9) indicate another that is composed of series 5 to 7 (corresponding to eigenvalue λ_{9}). In case 4, the participation indices PI(10) are equal to 1 (corresponding to the largest eigenvalue λ_{10}); this means that the cluster is composed of all time series and the contributions from each series are equal, i.e., the multivariate time series are completely correlated.

Figure 2 shows the relationship between the maximal (highest) eigenvalue and the number of synchronized elements of multivariate series. It is seen that the maximal eigenvalue increases linearly with the number of correlations (synchronizations) in the series when all the correlation strength is equal to 1. In brief, the closer the eigenvalue approaches the number of channels, the greater the degree of global synchronization between the multivariate series. Once this value is normalized it represents the global synchronization index; a value of 1 indicates that all the oscillators are fully synchronized and a value of 0 indicates that they behave as independent oscillators. Figure 2*B* shows the effect of channel number with respect to the synchronization index. The multivariate times series is composed of the correlated series generated by a Lorenz system (channels 1–4) and the noncorrelation series generated by a random system. The maximal eigenvalue is equal to the number of correlated series when their correlation coefficients approach 1 (Fig. 2*A*). With the addition of increasing numbers of random series, the synchronization index of the new multivariate time series decreases, as can be seen in Fig. 2*B*; however, the participation index remains the same (see Fig. 3). This means that caution needs to be exercised when comparing different recording configurations: electrodes placed outside the synchronized area will dilute the maximum value, but will have no impact on time-dependent changes, nor on the participation index, which identifies synchronized sites. Varying the correlation strength for case 4 in Fig. 1*A* shows that the maximal eigenvalue increases linearly with correlation strength (Fig. 2*C*).

Figure 3*A* shows a multivariate time series. Channels 1–4 are generated by a Lorenz system constructing a cluster; the other channels' time series are generated by a random system and are therefore not correlated. By changing the initial condition of the Lorenz system and randomly generating other channel time series, 100 multivariate time series may be obtained; the mean ± SD of the Syn_Index of the 100 multivariate time series is 0.245 ± 0.004. The participation indices of 100 multivariate time series are computed and plotted in Fig. 3*B*. The same procedure is applied to case 3 in Fig. 1, in which the multivariate time series is composed of two clusters: channels 1–4 and channels 5–7; other channel time series are generated by a random system, as can be seen in Fig. 3*C*. The mean ± SD of the Syn_Index of the 100 multivariate times series is 0.243 ± 0.004 for cluster 1 and 0.156 ± 0.003 for cluster 2, respectively. As can be seen in Fig. 3*D*, the participation index can clearly identify these two clusters in Fig. 3*C*. To perform a test for false positives, 100 multivariate time series are randomly generated; the mean ± SD of the Syn_Index is 0.071 ± 0.005, which is small in comparison with the results for the clusters in Fig. 3. These results show that this method can obtain a reliable synchronization index from multivariate time series and that elements of the cluster can be identified by the participation index.

To investigate the effect the noise on the above-cited method, first the multivariate time series in Fig. 3*A* are taken as pure signals, and white noise with different signal-to-noise ratios (0–10) are added to the pure signals. Then the Syn_Index and participation index are computed and the results are plotted in Fig. 4, *A* and *B*. As can be seen from Fig. 4*A*, the noise has a significant effect on the Syn_Index. Noise could be removed before the construction of correlation matrix, for instance using wavelet denoising (Li et al. 2007). In contrast, Fig. 4*B* shows that the participation index is relatively insensitive to noise.

### Analysis of multichannel epileptic field potential recordings: data recordings

Multiple field potentials were recorded, in vitro, from the CA1 region of the hippocampus. Transverse hippocampal slices (400 μm) were prepared from male Sprague–Dawley rats (150–300 g) and placed in an interface chamber. Spontaneous seizure-like events were induced by perfusion with low-calcium, high-potassium artificial cerebrospinal fluid consisting of (in mM) 125 NaCl, 26 NaHCO_{3}, 5 KCl, 0.2 CaCl_{2}, 1.0 MgCl_{2}, 1.25 NaH_{2}PO_{4}, and 10 glucose, bubbled with 95% O_{2}-5% CO_{2} mixture (Li et al. 2007). Extracellular field potentials were recorded using nine Pt/Ir wire electrodes (25 μm in diameter), positioned individually on the surface of stratum pyramidale along the long axis of CA1. The spacing of adjacent electrodes was about 200 μm. Signals were amplified using a Neuralynx headstage preamplifier (Neuralynx, Tucson, AZ), amplified (×500), low-pass filtered (3 kHz) with Neuralynx Lynx-8 amplifiers, and digitized with sampling frequency 5 kHz using 1401 Plus and Spike2 software (Cambridge Electronic Design, Cambridge, UK). The CMA and equal-time correlation methods were used to analyze the multichannel field potential recordings (band-pass filtered at 20–250 Hz). The purpose was *1*) to identify the number and size of synchronized clusters; *2*) to derive a synchronization index for each cluster; and *3*) to use the highest synchronization index as an indicator of global synchronization.

## RESULTS

### Comparison of three states

The correlations of nine channels are presented for three short multiple field potential (MFP) recordings; i.e., *M* = 9 and the length of of the analyzed epoch is 400 ms. The correlation matrix and its eigenvalues/surrogate eigenvalues are computed to show the different correlation (synchronization) structure of MFP recordings at the interictal and both early and late ictal states (Fig. 5). During the interictal state (Fig. 5*A*), the MFP reveals low-amplitude high-frequency activity (Fig. 5*Aa*), the correlation matrix shows no clusters (Fig. 5*Ab*), and none of the eigenvalues (Fig. 5*Ac*) exceeds the threshold based on distribution of surrogate eigenvalues obtained by 100 repeats using the AAFT method to destroy the phase relations of original MFP series.

During the early ictal phase (Fig. 5*B*) MFP recordings reveal larger activity on subsets of channels (Fig. 5*Ba*). The correlation (synchronization) matrix (Fig. 5*Bb*) shows one cluster, composed of channels 1–6. Comparison of the real and surrogate eigenvalues (Fig. 5*Bc*) confirms the existence of one cluster; the synchronization indices are 0, 0, 0, 0, 0, 0, 0, 0, 0.285. During the established ictal state the high-amplitude rhythmic activity is clearly visible on all channels (Fig. 5*Ca*). There is a single cluster composed of the channels 1–9 (Fig. 5*Cb*); the presence of one strong cluster is shown by the plot of eigenvalues (Fig. 5*Cc*). The synchronization indices are 0, 0, 0, 0, 0, 0, 0, 0, 0.836. The clustering strength is greater than that for the interictal and preictal states.

### From interictal to ictal state

Figure 6*A* plots MFP recording of nine channels during a 40-s period, starting during an interictal period and extending through a full ictal event. During epileptic seizures the activity of a large fraction of the neural network is strongly synchronized, and the MFP series displays an oscillation with large amplitude and regular frequency. When analyzing changes with time, the CMA analysis must be repeated. This means that an acceptable rate of false positives must be decided, and will determine the significance level used to calculate *K*, allowing for multiple comparisons as outlined earlier. Setting *K* = 3 with the nine channels used here results in a false positive rate of one per 100 epochs. In the present example, a window technique is applied to trace the synchronization of the MFP series over time. A moving window of 0.2 s with the overlap of 0.1 s is used. Figure 6*B* shows how the number of clusters varies between zero and two over the recorded period, with two clusters most prevalent during the early ictal stages. Plotting the two largest synchronization indices against time shows how the clusters appear and grow, and how their composition varies during the recording (Fig. 6, *C* and *E*): during the transition from the preictal to ictal period, the synchronization index corresponding to the largest eigenvalue (λ_{9}) increases gradually; the synchronization indices corresponding to the second largest eigenvalue (λ_{8}) grows during the early ictal period, and then declines as the seizure progresses. These results indicate that most channels appear synchronous toward the end of the recorded period. Figure 6, *D* and *F* shows the participation indices for the elements within the clusters (significant participation indices are represented by hot colors).

## DISCUSSION

It was previously shown that eigenvalue decomposition of a correlation matrix can quantify correlations in multivariate data sets and that this can be used as the basis for identifying clusters of synchronized oscillators (Müller et al. 2006). Herein, we have developed this method so that it can be applied to the analysis of neurophysiological data recorded simultaneously from multiple sites. Although synchronization of activity at two sites can be readily measured using existing techniques (Hurtado et al. 2004; Li et al. 2007), the present method enables analysis of the multiple site recordings that are frequently made, both experimentally (Feige et al. 2000; Fox et al. 2004; Taylor et al. 2005) and in clinical situations (D'Alessandro et al. 2005).

Selection of the most appropriate type of correlation used to construct the correlation matrix depends on the nature of the time series. When the time series is a smooth, continuous oscillation with very small (or zero) delays, equal-time correlation is simple and reliable; this is the correlation matrix we have used herein to illustrate the method. If the time series is interrupted by sudden bursts of irregular activity, equal-time correlation may still be used, but using the modification described by Müller et al. (2006); preprocessing (Govidan et al. 2005) can also be used to remove the effect the time delays between the time series. The matrix could, however, be composed from very different correlations, such as phase synchronization, linear cross-correlation, or event synchronization, depending on the nature of the activity being analyzed (Pereda 2005).

In principle, application of the method to neuronal data is very simple: a high eigenvalue indicates a high level of correlation and a low level indicates poor correlation. Because the maximum possible eigenvalue is equal to the number of recording channels, an eigenvalue approaching this figure indicates complete correlation of activity at all recording sites. When synchronization of neuronal activity is restricted to discrete areas, one or more relatively high eigenvalues are seen, corresponding to the number of areas showing synchronized activity. It is thus possible to identify clusters of local synchrony in which activity is synchronized and to determine whether that activity is independent of activity in other subpopulations within the recording area.

The synchronization index, derived from the largest eigenvalue, gives an indication of global synchronization throughout the population of neurons from which recordings have been made. It is established that epileptic seizures result from excessive synchronization of neuronal activity and we have therefore tested the procedure, using the low-calcium in vitro model of epilepsy. The results show that the technique can successfully identify synchronization of activity in small groups of neurons consisting of independent clusters of synchronous activity. Moreover, color coding of the participation index can be used to demonstrate how the distribution of clusters changes with time. The transition to an electrographic seizure is shown by an increase in the synchronization index as progressively more regions in the hippocampal slice are activated simultaneously.

Cognition and sensorimotor processing are frequently dependent on the synchronization of rhythmic activity in groups of neurons that may be anatomically circumscribed or may be widely separated through the neuraxis (Gray et al. 1989). Several methods have been successfully applied to analyze the interaction (synchronization) between multiple neural signals, including cluster analysis (e.g., Stuart et al. 2005), graph theoretic analysis (e.g., Stam et al. 2007), multichannel phase synchronization (Rudrauf et al. 2005), and mixture-of-Gaussians analysis (Matsumoto et al. 2005). However, each of these methods has its limitations and gives only partial information about the synchronization. Cluster analysis (including hierarchical cluster tree) often gives the information only about the connection architectures (the “tree”), but does not give information about the strength of the cluster. An interesting approach is conversion of the synchronization matrix into the binary graph and application of graph theoretic analysis. However, there is no unique way how to determine threshold and only one cluster can be described by this approach. Multichannel phase synchronization can reveal transient periods of ensemble synchronization in the time–frequency domain. However, it does not identify how synchronized clusters arise in an ensemble of oscillators, nor the degree of involvement of each oscillator in each cluster. The mixture-of-Gaussians analysis reveals the number of clusters and the elements of each cluster, but it does not give information about the strength of clustering. Like cluster analysis, it is not suited to track the interaction dynamics over time. Briefly, each of the existing methods has its limitations and disadvantages and they are not very well suited to track the collective dynamics of network during the time course, unlike the synchronization index and participation index in Fig. 6. We believe that the CMA analytical method presented here brings together all the relevant information needed to quantify synchronization between multiple time series.

## GRANTS

This work was supported by National Natural Science Foundation of China Grant 60575012 and Epilepsy Research Foundation (UK).

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society