Assessment of EEG synchronization based on state-space analysis

www.elsevier.com/locate/ynimg NeuroImage 25 (2005) 339 – 354

Assessment of EEG synchronization based on state-space analysis Cristian Carmeli,a Maria G. Knyazeva,b Giorgio M. Innocenti,c,d and Oscar De Feoa,* a

Laboratory of Nonlinear Systems, Swiss Federal Institute of Technology Lausanne, EPFL-IC-LANOS, Building EL E, Lausanne CH-1015 Switzerland Centre Hospitalier Universitaire Vaudois, Lausanne, Switzerland c Department of Neuroscience, Karolinska Institutet, Sweden d Department of Child Psychiatry, Centre Hospitalier Universitaire Vaudois, Lausanne, Switzerland b

Received 17 January 2004; revised 6 October 2004; accepted 30 November 2004

Cortical computation involves the formation of cooperative neuronal assemblies characterized by synchronous oscillatory activity. A traditional method for the identification of synchronous neuronal assemblies has been the coherence analysis of EEG signals. Here, we suggest a new method called S estimator, whereby cortical synchrony is defined from the embedding dimension in a state-space. We first validated the method on clusters of chaotic coupled oscillators and compared its performance to that of other methods for assessing synchronization. Then nine adult subjects were studied with high-density EEG recordings, while they viewed in the two hemifields (hence with separate hemispheres) identical sinusoidal gratings either arranged collinearly and moving together, or orthogonally oriented and moving at 908. The estimated synchronization increased with the collinear gratings over a cluster of occipital electrodes spanning both hemispheres, whereas over temporo-parietal regions of both hemispheres, it decreased with the same stimulus and it increased with the orthogonal gratings. Separate calculations for different EEG frequencies showed that the occipital clusters involved synchronization in the b band and the temporal clusters in the a band. The g band appeared to be insensitive to stimulus diversity. Different stimulus configurations, therefore, appear to cause a complex rearrangement of synchronous neuronal assemblies distributed over the cortex, in particular over the visual cortex. D 2004 Elsevier Inc. All rights reserved. Keywords: S estimator; Cortical synchrony; State-space

Introduction The operations performed by the cerebral cortex require the cooperative activity of neurons within distributed assemblies. In human studies, the assemblies have been identified with fMRI, which provides a detailed static image but poor temporal resolution, with EEG techniques, which provide a dynamic image but poor spatial resolution, and more recently, with a combination

* Corresponding author. Fax: +41 21 69 36700. E-mail address: [email protected] (O. De Feo). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2004.11.049

of those two methods (see, for instance, Goldman et al., 2000; Knyazeva et al., 2002, 2003; Lemieux et al., 2001; Mulert et al., 2002). The use of EEG in the studies of neuronal assemblies has been guided by several important concepts. The first one is that at least certain types of assemblies are characterized by the synchronous activity of their constituent neurons. The second one is that the different EEG frequency components also reveal synchronies relating to different perceptual, motor or cognitive states (Anokhin et al., 1999; Basar et al., 2001; Burgess and Gruzelier, 1997; Eckhorn et al., 1988; Engel et al., 1991; Klimesch, 1999; Klimesch et al., 1994; Miltner et al., 1999; Rodriguez et al., 1999; Tallon-Baudry et al., 2001; Tallon et al., 1995). These two concepts are deeply rooted in the tradition of the EEG community as well as in the recent findings of animal and human neurophysiology. They have been supported by a lot of studies showing, for example, increased power in the b–c range of EEG signals during perceptual or visuomotor tasks (Aoki et al., 1999; Keil et al., 1999; Muller et al., 1996). The most common way of analyzing the cooperative, synchrony-defined cortical neuronal assemblies has been the study of EEG coherence. This method identifies the synchrony of neuronal assemblies as a function of the correlation of EEG frequency components. Neurons in an assembly are presumably located in some proximity to the recording electrodes and exhibit oscillatory activity with common spectral properties. The typical finding is that in a given perceptual, cognitive, or motor task, EEG coherence increases (or decreases) in a certain band of EEG frequency spectrum (Ford et al., 1986; Kiper et al., 1999; Knyazeva et al., 1999; Mima et al., 2001). It is usually unclear if the functional connectivity revealed by the increased EEG coherence signifies anatomical connectivity between the neurons of the assembly (Friston, 1994; Friston et al., 1993; Horwitz, 2003), although this assumption has been verified in a very restricted number of experimental and/or clinical conditions (Koeda et al., 1995; Kuks et al., 1987; Lopes Da Silva et al., 1980; Montplaisir et al., 1990; Pinkofsky et al., 1997). The definition of synchrony used in the neurophysiological studies is narrower than that considered in system theory and, in particular, in the study of complex systems (Boccaletti et al., 2002;

340

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

Brown and Kocarev, 2000). Indeed, while the definition used in neurophysiology assumes two or many subsystems sharing specific common frequencies, in system theory bsynchronizationQ refers to a process whereby two or many subsystems adjust some of their timevarying properties to a common behavior due to coupling or common external forcing. Hence, the general definition of synchronization does not imply a relationship between systems at specific frequencies but at all of them. Clearly, the two above definitions coincide in the case of perfectly periodic systems, though coherence at specific frequencies is not sufficient to fully define synchrony in a system theory sense. A figurative example will help in understanding the occurrence of non-periodically synchronized systems. Let us consider a symphonic orchestra. During a performance, each player is continuously synchronized with the others by means of either the action of the Conductor (common forcing) or by listening to each other (mutual coupling). However, the symphony is neither strictly nor approximately periodic. Although the methods focusing on the analysis of synchronization phenomena restricted to particular frequency bands have been the norm in EEG studies, work grounded on a system theory definition of synchronization has begun to appear (Quiroga et al., 2002; Schiff et al., 1996; Tass et al., 2003). Unfortunately, the assessment of synchronization phenomena can be difficult if one adheres too strictly to the mathematical definition; consequently, most of the success of these new studies has been in cases where the synchronization manifests strongly, for instance, in sleep analysis (Fell et al., 2003), epileptic seizures (Mormann et al., 2000), and Parkinsonian subjects (Tass et al., 1998). On the other hand, quantifiers of synchronization resorting to thermodynamics as well as to information theory can be defined, such as entropy and mutual information (Inouye et al., 1991; Quiroga et al., 2000; Sabesan et al., 2003). Being robust, these quantifiers can be applied to more general EEG contexts (Stam et al., 2003). Along this line, in this paper, we have quantified the synchronization of EEG activity induced by visual stimuli by measuring the contraction of the embedding dimension of the statespace, which is a theoretical consequence of synchronization (Boccaletti et al., 2002). The stimuli were initially designed, and successfully used, to identify the synchronization of activity between the hemispheres (Knyazeva et al., 1999). Therefore, the results we obtained could be compared with those of the more traditional EEG coherence method, employed in that study, and are in general consistent with the latter. However, the complexity of the changes induced by those simple stimuli seems to go well beyond what we had expected from previous work.

Materials and methods The experimental conditions are identical to those used in (Knyazeva et al., 1999, 2002, 2003). The data analysis differs since in this paper we introduce a method based on the contraction of the embedding dimension (Boccaletti et al., 2002), in order to estimate the degree of regional synchronization of cortical activity. Experimentalia We have used the EEG recordings from nine normal righthanded adults (7 women, 2 men; mean age 30 years, range of 27–

47 years), who participated in experiments aimed at mapping neural assemblies induced by visual stimulation (Knyazeva et al., 2002, 2003). This had the advantage of testing our method with an experimental paradigm known to result in reproducible synchronization dynamics, measured with spectral coherence analysis. All the subjects, CD, EF, FG, KL, MR, NM, PH, UH, and VO signed written informed consent and all the procedures conformed to the Declaration of Helsinki (1964) by the World Medical Association concerning human experimentation and were approved by the local ethical committee of Lausanne University. During the EEG recording session, the subjects viewed the visual stimuli shown in Fig. 1. The stimuli were generated with a PC and presented on the computer display with a refresh rate of 75 Hz located at a distance of 57 cm. The subjects were instructed to fix a point in the center of the screen. The stimuli were black and white bilateral iso-oriented or orthogonally oriented sinegratings centered on the fixation point. The iso-oriented gratings (IG) consisted of two identical patches of collinear, downwardsdrifting horizontal gratings on both sides of the fixation point. The orthogonally-oriented gratings (OG) consisted of a patch of horizontal downwards-drifting grating on one side and a patch of vertical rightwards-drifting grating on the other. All the gratings had a spatial frequency of 0.5 cycles/degree, a contrast of 70%; the unilateral patches measured 13.58 (width) by 248 (height). They drifted with a temporal frequency of 2 Hz. A uniform gray screen of the same space-averaged luminance as the stimuli (32 cd/m2) with a fixation point in the center served as a background stimulus. The vertical and horizontal gratings of the OG stimulus appeared in the left or right hemifield at random. The type of stimulus (OG, IG, and background), the stimulus exposure (2.2–2.6 s), and the inter-stimulus intervals (1.8–2.2 s) were also randomized. The EEGs were recorded with a 128-channel Geodesics Sensor Netn. All electrode impedances were kept under 50 KV. The recordings were made with vertex reference using a low-pass filter set to 100 Hz. The signals were digitized with a 12-bit analog-todigital converter at a rate of 500 samples/s. They were further filtered (FIR, band-pass of 3–70 Hz, notch of 50 Hz), rereferenced against common average reference, and segmented into non-overlapping 1 s epochs using NS2/NS3 (Electrical Geodesics, Inc.) software. With adequate sampling of the head surface provided by our 128-channel EEG system, the average-reference EEG signals provide a reasonable approximation of referenceindependent estimates of the EEG sources (Srinivasan, 1999, 2003). Artifacts were edited off-line, first automatically, based on an absolute voltage threshold (100 AV) and on a transition threshold (50 AV), and then through visual inspection, which allowed the identification of channels and/or epochs with moderate muscle artifacts not reaching the threshold values. Finally, in order to assess only the steady-state, the first 200 ms following stimulus onset were removed, thus excluding stimulus-onset artifacts, response-onset transients, and stimuluslocked synchronization. Furthermore, in order to estimate synchronization independently from the absolute energy of the signals, each signal was mean de-trended and scaled to unit time variance. We collected on average 85 artifact-free epochs for each stimulation condition and for each individual, with a minimum of 78 and a maximum of 97.

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

341

Fig. 1. The four different types of visual stimuli.

The S estimator The term bsynchronizationQ colloquially refers to time-correlated behavior between different processes. However, a standard, unique, rigorous definition of synchronization does not exist and particular definitions are adopted according to the application or the investigation considered (Brown and Kocarev, 2000). As mentioned in Introduction, we define as synchronization a process whereby two (or many) dynamical subsystems adjust some of their time-varying properties to a common behavior due to coupling or common external forcing. More precisely, considering a large stationary, deterministic, finite dimensional dynamical system, which can be divided into two n and m dimensional subsystems, respectively 8 dx > > < dt ¼ f 1 ðx;y;t Þ x a Rn ; y a Rm ð1Þ > dy > : ¼ f 2 ðx;y;tÞ dt we say that two sub-trajectories x(t) and y(t) of the whole system are synchronized with respect to the properties (time dependent measures) gx and gy  gx : Rn  R Y Rk k V minðm; nÞ ð2Þ gy : Rm  R Y Rk if there is a time independent mapping h: Rk  Rk Y Rk such that  h i   ð3Þ h gx ðxÞ; gy ðyÞ  ¼ 0 where jj!jj is any norm in Rk . As explained by Brown and Kocarev (2000), this unifying definition covers most kinds of usually considered synchronizations, such as identical synchronization, phase locking synchronization, and phase synchronization. In brief, condition (3) requires a property (gx) of the trajectory x(t) to be in a fixed relation (h) with another property (gy) of the trajectory y(t). As an example, let us consider two periodic trajectories (not necessarily equal) with the same frequency: they are synchronized

according to condition (3) if we assume gx,y to be the frequency and h = gx  gy. Actually, condition (3) implies that synchronized sub-trajectories lay on an r-dimensional manifold, where r depends on h and is 1 V r V k. Consequently, the dimensionality of synchronized dynamics (n + m  r) results to be smaller than that of the generic asynchronous dynamics (n + m) in the whole system (Brown and Kocarev, 2000; Boccaletti et al., 2002). In other words, let us consider two independent dynamical systems X and Y, whose state-spaces are of dimension n and m, respectively. The compound system obtained considering the two systems together has dimension n + m. However, the embedding dimension of the dynamics of the compound system turns out to be smaller than n + m when the two subsystems are synchronized. For instance, two identically synchronized 2-dimensional periodic oscillators have dimension of the whole state-space 4, but the embedding dimension of the synchronized dynamics is only 2. Exploiting what has just been discussed, we propose to quantify the amount of synchronization within a data set by comparing the actual dimensionality of the set with the expected full dimensionality of the asynchronous set. We propose to perform this comparison by considering the embedding technique based on Principal Components Analysis (PCA) (Broomhead and King, 1986; Kantz and Schreiber, 1997)1. The proposal of a new measure based on the PCA embedding technique, rather than exploiting one of the synchronization estimators existing in literature, follows from the multivariate nature of EEG recordings. As it is clear from its definition, synchronization is a bivariate concept. Consequently, the direct synchronization estimators proposed in literature (e.g., Quiroga et al., 2002) are intrinsically bivariate and hard to extend to multivariate cases. However, synchronization should be often 1 An alternative justification of this synchronization descriptor has been proposed by (Wackermann, 1996) drawing inspiration from thermodynamics.

342

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

being proportional to the embedding dimension of the dynamical phenomenon. In fact, the more disperse the eigenspectrum is, the more numerous the significant PCA components are (which means higher embedding dimension), and the higher the entropy of the eigenspectrum. On the contrary, the more concentrated the eigenspectrum is, the fewer the significant PCA components are (which means lower embedding dimension), and the lower the entropy of the eigenspectrum. Though this measure is information theory inspired, it is not rigorously entropy, as the normalized eigenspectrum is not a probability distribution. Let us consider K-variate non-delay-embedded EEG recordings; the same holds, mutatis mutandis, for the case of delayembedded signals. Each measurement is composed of voltage vectors u i that can be collected in a matrix U such that 2 3 u0 6 u1 7 K 7 ð5Þ U¼6 4 v 5ui a R uN 1

estimated among multivariate time series, which is the case in the analysis of EEG recordings. The use of PCA, which is a multivariate technique, provides an intrinsically multivariate method. Indeed, PCA explains the variance–covariance structure of multivariate data through a few linear combinations of the original variables. A given multivariate time series with k components can be transformed into the so-called population principal components (PC) by a linear transformation projecting the original time series into the eigenbase of the covariance matrix (correlation matrix, if the data have been power normalized) of the time series itself (Joliffe, 2002). In this coordinate system, the relative importance of each principal component, in justifying the variance of the original time series, is given by the normalized eigenvalue associated to the eigenvector corresponding to the considered component. In this way, the PCA is usually adopted for data reduction. In fact, although the original data set contains k variables, often much of the variability can be accounted for by a smaller number, m b k, of principal components. Accordingly, PCA is often used to determine the embedding dimension and the embedded trajectory of a dynamical phenomenon (Broomhead and King, 1986; Kantz and Schreiber, 1997). This is usually carried out by means of the delay-embedded (DE) coordinates. Namely, a given univariate time series y(t) can be transformed into a multivariate one through the delay embedding

where K is the number of considered electrodes and N is the number of measurements for each epoch. The corresponding K K correlation matrix is given by 1 1 NX C¼ un uTn ð6Þ N n¼0

yðt Þ Y HðtÞ ¼ ½ yðt Þ yðt  sÞ : : : yðt  KsÞ

If k 1,. . ., k K are its eigenvalues, and

ð4Þ

where the time delay s and the window size K may be computed through the autocorrelation or the self-mutual information (Kantz and Schreiber, 1997). Thus, as Q(t) has the form of a multivariate time series, PCA can be applied: the embedding dimension is given by the number of principal components of the DE coordinates necessary to justify the variance of the dynamical phenomenon itself and, on the other hand, the embedded trajectory is obtained by projecting the DE coordinates into the principal components. This technique is trivially extended to the case of multivariate time series by applying the delay embedding to each component. The delay embedding is formally necessary to reconstruct the state-space embedded trajectory of the dynamical phenomena under observation from the measured scalar time series (Kantz and Schreiber, 1997). However, for multivariate time series, this is computationally and memory expensive and it may be avoided when considering highdensity surface sampling EEG recordings. Indeed, for such recordings, closely positioned electrodes probably bring information about the same oscillator (state-space) providing an already embedded trajectory of the dynamical phenomenon, as further corroborated by a complex numerical experiment reported in Assessing the topographical sensitivity and robustness of S. Hence, in the remainder of this paper, we will consider k-variate non-delay-embedded EEG recordings; although it is to be noticed that the method can be applied in a similar fashion to delay-embedded data. Considering the case of non-delay-embedded multivariate signals, it is then possible to define a measure of the amount of synchronization, for a region of the cortex, relying solely on the eigenspectrum of the correlation2 matrix of a group of electrodes. Our proposal is to measure the amount of synchronization over a region of the cortex by an information theory inspired measure defined as the complement of the entropy of the normalized eigenvalues of the corresponding correlation matrix; such quantity 2 The covariance matrix is in reality the correlation one, as the signals have been normalized in variance as explained in Experimentalia.

ki V ¼

ki K P

¼ kj

ki trðCÞ

ð7Þ

j¼1

the corresponding normalized eigenvalues, the quantity K P

S ¼1þ

ki Vlogðki VÞ

i¼1

logð K Þ

ð8Þ

is a measure inversely proportional to the embedding dimension of the observed dynamical phenomenon, thus proportional to the amount of synchronization. To understand how this measure works, let us consider first the case with a single field generator for the entire skull. Under this condition, the correlation matrix C degenerates to the case with k 1V = 1 while all other eigenvalues are zero, so that S = 1. On the other extreme, if there are K uncorrelated generators under each electrode site, the correlation matrix C is diagonal (it is in fact an identity matrix because of signals’ normalization) and all the normalized eigenvalues are expected to be kVi = 1/K, so that S = 0. Note that in this sense, in the case of non-delay-embedded signals, the S estimator can be interpreted as a deviation from mutual orthogonality (lack of correlation) among the K EEG recordings. In fact S is equal to 0 if there are K uncorrelated (orthornormal because of power normalization) signals, and S is equal to 1 if they are perfectly correlated. However, this interpretation does not hold in the case of delay-embedded signals. Furthermore, it should be noted that the definition of S is always independent from both the total power and the time dimension; hence, S is independent from any time-based parameter such as frequency. Furthermore, for the case of non-delayembedded signals, it inherits from the eigenvalues of the PCA the rotational invariance, that is, S does not change with the coordinate system: time, frequency, or any other orthogonal base (Joliffe, 2002). These two extreme cases explain how S works as an estimator of spatial synchronization over the array of K electrodes. However,

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

they do not illustrate the relationship between S and the classical measure used in electroencephalography, the coherence. This relationship can be demonstrated for bivariate signals. In this simple case, the correlation matrix C is given by   1 q C¼ ð9Þ q 1 where q is the correlation between the two signals. Consequently, the normalized eigenvalues are k 1,2 V = ((1 F q) / 2) and S is given by S¼

ð1  qÞlogð1  qÞ þ ð1 þ qÞlogð1 þ qÞ 2logð2Þ

ð10Þ

which is a nonlinear monotonically increasing convex function of the correlation q. By the Parseval’s theorem, the correlation q is the integral over all the frequencies of the coherence, showing that in this sense the S estimator is spectrum wide and not frequency specific. It is clear that in the case of multivariate or delay-embedded (state-space extended) signals, the relationship is not so trivial and an analytical analysis cannot be carried out easily3. In Numerical assessment of the S estimator, the relationship between the estimator S and the coherence, as well other measures of synchronization, is illustrated together with a robustness assessment. Application to high surface sampling EEG Previous studies based on coherence analyses (Knyazeva et al., 1999, 2002) have shown that, for the stimuli used, an interhemispheric synchronization, probably implemented via callosal connections, can be detected at occipital electrodes close to the midline (e.g., electrode pairs 71–84 and 70–90 in Fig. 2). In the present study, the phenomenon of interhemispheric synchronization was also of interest and it was estimated that the adequate bihemispheric territory was covered within a distance of one or two electrodes across the midline. This corresponds to a distance of 6–9 cm and is reasonably well spanned by the second order neighborhood of each electrode (cf. Fig. 2). On the contrary, the first order neighborhood would be too small to credibly span across the two hemispheres, as shown in Fig. 2, while the third order one would be too large, including about half of the electrodes at once. Therefore, we considered the second order neighborhood assuming that contiguous regions of synchronization (positive DS) crossing the midline could be interpreted as highlighting interhemispheric synchronization. However, it is clear that it would have been possible to assess distant synchronization phenomena by computing S for sets of far electrodes on the left and right hemispheres. To summarize, the close interhemispheric synchronization was assessed by computing, for each of the nine subjects mentioned in Experimentalia, the topographical arrangement of the synchronization activity as follows. First, for each electrode, the synchronization of the surrounding region was computed as the S value over its

3

Actually, for the case of multivariate (non-delay-embedded) data, it is possible to show that S is monotonically increasing with respect to all the correlation terms. The proof is based on the diagonalization of a symmetric matrix by means of Jacobi rotations (Press et al., 1992), and it is available from the authors upon request.

343

second nearest topographical neighbors (Fig. 2) and an instantaneous topographical map of synchronization values (S) was obtained for each epoch. The instantaneous values of epochs corresponding to identical stimulus condition (BGR, IG, and OG, cf. Experimentalia) were averaged obtaining the statistics, mean vector and covariance matrix, of the synchronization topography corresponding to each one of the three visual stimuli. Then, to highlight the reorganization of the topographical activity in the presence of different stimuli with respect to the background activity, we subtracted the average value of the S estimator in the background from the average value in the two actively stimulated situations. In what follows, we refer to this difference as DS = S IG,OG – S BG , where its positive values denote increased synchronization with respect to the background, and negative values indicate decreased synchronization. Therefore, extracting the regions corresponding to the maxima of this differential synchronization estimator, we have highlighted the clusters of electrodes with synchronized behavior. Statistical assessment of the results As explained in the previous section, in order to evaluate the topographical arrangement of stimulus-induced synchronization, we consider the variation DS of the mean values of the S estimator in the presence of the two grating conditions with respect to the values in the background. Hence, it is necessary to evaluate the statistical significance of DS. The values of S over the electrode array constitute a multivariate population. Usually, for multivariate data the significance of the difference between the means of two populations is addressed with the one-way MANOVA (Hotelling’s T 2) test (Johnson and Wichern, 2002). However, when the hypothesis of a normal distribution of the data is not tenable, as in our case (S a [0,1]), one may resort to the permutation version of the Hotelling’s T 2 (Higgins, 2004). This version consists in a T 2 test, where the reference statistic distribution F is substituted by the empirical distribution of the randomly permuted data. Accordingly, we have verified the statistical significance of DS, for individual cases as well as for the population, by considering the permutation version of the Hotelling’s T 2 without the assumption that the two covariance matrices are equal (Krishnamoorthy and Jianqi, 2004). For the individuals, the statistics were restricted to a group of artifact-free electrodes within a region of interest, and for the population, to the electrodes, which were artifact-free for all the subjects (cf. Results).

Numerical assessment of the S estimator Since it is not easy to assess analytically the properties and performances of the S estimator, we have proceeded with a numerical assessment. In particular, we wanted to evaluate three aspects of the S estimator: (i) that it measures indeed synchronization, and to compare it with other classical as well as innovative measures of synchronization; (ii) its robustness with respect to observation noise and length; (iii) its topographical sensitivity and robustness. The first two aspects were addressed considering bivariate-simulated signals, and the third by considering multivariate signals caricaturing the experimental EEG conditions.

344

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

Fig. 2. Example of the spatial localization of the synchronization estimator S. As an example, first nearest neighbors (in gray) are shown for electrode 76 (in green), first and second neighbors (in gray) for electrode 28 (in green). The latter has been used in the computations reported here.

Assessing S as a measure of synchronization and its robustness In order to determine how the S estimator depends upon synchronization, we have considered a classical synchronization paradigm (Boccaletti et al., 2002; Quiroga et al., 2002): two different nonlinearly coupled chaotic oscillators, that is, a Rfssler oscillator (Ro¨ssler, 1976), nonlinearly driving a Lorenz one (Lorenz, 1963), as given by 8 dx1 > > ¼  T ðx1 þ x3 Þ > > dt > > > >dx2 > > ¼ T ðx1 þ ax2 Þ > > dt > > > > > dx3 > > > < dt ¼ T ½bx1 þ x3 ðx1  cÞ ð11Þ dy1 > > ¼ rðy2  y1 Þ > > > dt > > > > dy2 > > ¼ qy1  y2  y1 y3 þ C ðx2  y2 Þ3 > > > dt > > > > dy3 > > ¼ y1 y2  by2 : dt where x i are the state variables of the Rfssler oscillator, y i those of the Lorenz one; a, b, c, b, r, and q are parameters fixed at the standard values, that is, a = 0.2, b = 0.2, c = 5.7, b = (8/3), r = 10, q = 28, the time scale T = 6 adapts the bspeedQ of the Rfssler to the

one of Lorenz, and the parameter C is varied in the range [0,1] so that to change the degree of synchronization between the two systems. For every considered value of C, the differential Eq. (11) were iterated, starting from random initial conditions, using the Heun algorithm (Quarteroni et al., 2000) with Dt = 0.005, which were checked to yield numerically stable results. In order to eliminate transients, the first 10000 iterations were discarded, then after performing a down-sampling operation at DT = 0.02, time series of length N = 500, the length of one epoch in our EEG data, were collected. To ensure repeatability and stability of results, we collected a total of 80 epochs from different initial conditions. We measured the coupled variables x 2 and y 2, and we added zero mean white Gaussian observation noise of different intensities, that is, 40, 20 and 12 dB Signal to Noise Ratio (SNR). Afterwards, we assessed the relationship between synchronization and the estimator S, as well as the robustness of S with respect to measurement noise and data length. This was done for each of the values of C and the observation noise intensities, by comparing four different measures of synchronization computed over 20, 40, or all the 80 epochs. The measures considered were (Quiroga et al., 2002): coherence, mutual information, the H estimator and our proposed S estimator. Excluding coherence, for which this is not possible, all the synchronization measures were computed considering the delay-embedded multivariate signals. More precisely, we delay-embedded both the measured time series in a 3dimensional space using embedding delays of s x = 14 samples

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

for x 2 and s y = 12 samples for y 2; these delays were obtained through the minima of the self-information (Kantz and Schreiber, 1997). Hence, we obtained the multivariate time series ui ¼ ½x2 ðiDT Þ x2 ðði  700ÞDT Þ

x2 ðði  1400ÞDT Þ

y2 ðiDT Þ y2 ðði  600ÞDT Þ y2 ðði  1200ÞDT Þ

ð12Þ

for which we computed S with the technique described in The S estimator. The results are summarized in Fig. 3, which compares the dependence of the four synchronization estimates upon the coupling parameter at the different noise intensities and the number of epochs considered. S scales similarly to the other estimators with respect to the coupling parameter and it similarly falls at C g 0.1 because of a bifurcation of the driven oscillator. For all the estimators, the computed values decrease when noise is increased; nevertheless, S appears more robust since it decreases less than the others. Furthermore, as expected, the reliability of S (cf. the error bars) increases with data length, similarly to the other estimators. As far as coherence is concerned, the high values at very low coupling suggest that coherence might reflect similar spectral properties of the two oscillators rather than real coupling (connectivity). Finally, the reader should not be surprised by the fact that S does not scale from 0 to 1 as C varies from 0 to 1. Indeed, the

345

whole state-space is 6-dimensional and the two isolated strange attractors have fractal dimension g2.01 and g2.06 for Rfssler and Lorenz, respectively. Hence, for C = 0, the dynamical phenomenon observed has fractal dimension g2.01 + 2.06 = 4.07 which may be contained in a 5-dimensional space, consequently S = 1  5/6 g0.17 p 0, which is indeed similar to the measured minimum. On the other hand, when C = 1, if the two systems were identically synchronized, S would be S = 1  q2.01a / 6 = 0.5 p 1, as only the Rfssler oscillator would determine the dimensionality of the phenomenon. However, for C = 1, the two systems are strongly synchronized albeit not identically as shown by the value of S, which is slightly lower than predicted. Assessing the topographical sensitivity and robustness of S If we should consider just bivariate data, there would be no need to introduce a new estimator of synchronization, since any of those considered above would suffice. However, the estimator S is interesting for its extendibility to multivariate data. To assess the topographical sensitivity and robustness of the S estimator, we resorted to simulated data approaching the situation described in Experimentalia. We considered the simulation of a

Fig. 3. Numerical assessment of the S estimator as a measure of synchronization. Dependence of the synchronization between the subsystems x and y of system (11), upon the coupling parameter C, estimated with different techniques: S estimator (red); mutual information (green); H estimator (cyan); and maximum of spectral coherence (blue). Robustness of the synchronization measures with respect to observation noise, considering 40, 20, and 12 dB Signal to Noise Ratio, and data lengths, considering 20, 40, and 80 epochs of 500 samples each. The dots and errors bars illustrate the mean value and standard deviation, respectively.

346

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

network of coupled 128 chaotic Colpitts oscillators. The Colpitts is a chaotic nonsymmetrical oscillator which generates irregular sinelike signals such as the EEG (De Feo et al., 2000); hence, we preferred it to the Lorenz and Rfssler oscillators considered above. The dynamics of the oscillator at site i is given by 8 ðiÞ h  ðiÞ i  dx1 ð iÞ g > x > > dt ¼ Qð1  k Þ a e 2  1 þ x3 > > > 128 h i X    ðiÞ  < ðiÞ dx2 ð iÞ ð jÞ ð iÞ g x2 þ ¼ ð 1  a Þ e  1 þ x C x  x i;j 3 2 2 dt Qk > > j¼1 > >   ðiÞ > > : dx3 ¼  Qk ð1  k Þ xðiÞ þ xðiÞ  1 xðiÞ 1 2 dt g Q 3 ð13Þ where g, Q, a, and k are parameters. In the simulations, k = 0.5 for all the oscillators, while g, Q, and a have been chosen different for each oscillator, with values in the intervals [4.006, 4.428],

[1.342,1.483], and [0.949, 0.999], respectively. Consequently, the network was composed of 128 non-identical chaotic oscillators (De Feo et al., 2000). The rightmost sum term on the right hand side of x 2(i) describes the coupling from the other oscillators j to the oscillator i. We considered three different coupling schemes assuming a very simplified hypothesis about the three stimuli situations BGR, OG and IG. Respectively: (i) all oscillators being isolated (BGR), that is, C i,j = 0 8 i, j; (ii) two isolated ring-clusters (OG) of three coupled oscillators each: a cluster of electrode 65, 66, and 71, and another (symmetric) one of electrodes 84, 85, and 91, that is, C i, j = 0 8 i, j except C 66,65, C 71,66, C 65,71, C 85,91, C 84,85, and C 91,84 which are nonzero and positive; and, (iii) the two clusters mentioned above unidirectionally coupled (IG) setting C 84,71 to nonzero and positive. In the simulations, the nonzero, C i, j were randomly chosen (once for all) in the range [1,3], which guarantees strong coupling but not perfect synchronization.

Fig. 4. Numerical assessment of the topographical sensitivity of the S estimator. Topographical maps of the variations of the mean synchronization (DS, see text) of the simulated EEG activities OG vs. BGR. Robustness of the synchronization topography with respect to observation noise, considering 10, 4, and 0 dB Signal to Noise Ratio, and data lengths, considering 20, 40, and 80 epochs of 500 samples each. White-bordered dark regions highlight the clusters of electrodes with the highest DS values. Five white contour lines emphasize the rate of transition from synchronized to unsynchronized clusters. The contour values are equally spaced between DS = 0 (no augmentation of synchronous activity with respect to the background) and the maximal value of DS. The shade of colors is obtained by a tri-linear interpolation from the three nearest electrodes. Gray and red spots are the electrodes positions; the red spots show the position of the two synchronized clusters considered in the simulation (see text).

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

Similarly to what is described above, the three networks were simulated with the Heun method (Dt = 0.025) starting every time from random initial conditions. Then, the transients were dropped out and the data down-sampled (DT = 0.075) collecting 80 epochs of 500 points each. Again, we measured the coupled variables x 2(i) from all sites, and we added zero mean white Gaussian observation noise of different intensities, that is, 10, 4, and 0 dB SNR. Afterwards, we assessed the topographical sensitivity and the robustness of the S estimator with respect to measurement noise and data length. This was done for each of the values of the observation noise intensities analyzing the data as if they were real EEG according to the procedure described in Application to high surface sampling EEG, considering at once only 20, 40, or all the 80 epochs. The results are summarized in Figs. 4 and 5, which show the topography of DS values for the IG vs. BGR and OG vs. BGR arrangements, respectively, at the different noise intensities and number of epochs. In both figures, the white-bordered dark regions highlight the territory with the highest DS values. In order to emphasize the steepness of the transition from synchronized to

347

unsynchronized clusters, five white contour lines are drawn. The contour values are equally spaced between DS = 0 and the maximal value of DS. Referring to the best case (80 epochs and 10 dB SNR), it can be noticed that the synchronized clusters are correctly identified by the topography of DS without delay-embedding, even if in this complex case each signal comes from a different oscillator. Furthermore, the second-order neighborhood cluster analysis, as suggested in Application to high surface sampling EEG, is indeed adequate to observe close interhemispheric synchronizations, possibly caused by the IG stimuli, as shown by the wide crosshemispheric synchronized regions in Fig. 5. On the contrary, when the clusters are uncoupled, the topography of DS also shows isolated regions of synchronization, cf. Fig. 4. As the SNR decreases, the synchronized clusters visually smear out at the chosen equal plotting scales. Nevertheless, the synchronized clusters, as well as the difference between the two topographical arrangements, can be verified to be statistically significant down to 4 dB with 40 epochs and to 0 dB with 80 epochs, as it was the case for the real EEG.

Fig. 5. Numerical assessment of the topographical sensitivity of the S estimator. Topographical maps of the variations of the mean synchronization (DS, see text) of the simulated EEG activities IG vs. BGR. Robustness of the synchronization topography with respect to observation noise, considering 10, 4, and 0 dB Signal to Noise Ratio, and data lengths, considering 20, 40, and 80 epochs of 500 samples each.

348

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

Results The EEG data from nine subjects have been analyzed with the synchronization estimator S. This estimator highlighted regions whose activity was synchronized by the visual stimuli. Further, we broke down the analysis of the S estimator to three of the traditional frequency bands of EEG analysis, the a (7–13 Hz), b (13–30 Hz), and c (30–70 Hz) bands. Synchronization topography The visual stimuli used in this study were chosen primarily because it was previously found that they could reliably modify interhemispheric spectral coherence in animals and men (Knyazeva et al., 2002, 2003). Stimulus IG conforms to the Gestalt principle of perceptual grouping by collinearity and common fate and is likely to generate cooperative activity in the two hemispheres. Stimulus OG, instead, presumably places the visual area of the two hemispheres in a condition of perceptual rivalry. Consequently, according to the synchronization-based solution to the binding problem (Gray, 1999), high S

values are predicted in areas of the brain involved in processing the IG stimulus, while the rest of the activity, remaining unorganized, should result in low S values. Indeed, as explained above, a high value of S reflects a low amount of independent activity sources in the region processing the coherent stimulus. The average topographical arrangement of DS values with the IG and OG stimuli are shown in Figs. 6 and 7. In both figures, the white-bordered dark regions denote the clusters of electrodes (and indirectly the underlying cortical territory) with the highest DS values. As for the simulations (cf. Assessing the topographical sensitivity and robustness of S), in order to render the steepness of transition from synchronized to unsynchronized clusters, five white contour lines are superimposed to the figures. The contour values are equally spaced between DS = 0 (no augmentation of synchronous activity with respect to the background) and the maximal value of DS. Among the nine subjects, the region of synchronization induced by the OG stimulus (Fig. 6) is more often lateralized, while the one induced by the IG stimulus (Fig. 7) is more often centered at or near the electrodes close to the occipital lobes and

Fig. 6. Topographical maps of the variations of the mean synchronization (DS, see text) of EEG activity induced by the orthogonally oriented grating (OG) with respect to the background. White-bordered dark regions highlight the clusters of electrodes (and indirectly the underlying cortical territory) with the highest DS values. Gray and red spots are the electrodes positions, the red ones being excluded from the computation due to artifacts. CD, EF, FG, KL, MR, NM, PH, UH, and VO denote the subjects. Other conventions as in Fig. 5.

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

349

Fig. 7. Topographical maps of the variations of the mean synchronization (DS, see text) of EEG activity induced by the iso-oriented grating (IG) with respect to the background. Other conventions as in Fig. 6.

spans the two hemispheres. This is further highlighted in Fig. 8, which shows the response to the IG vs. the OG stimulus. Nevertheless, in this and in previous studies using the same stimuli and EEG coherence assessments (Knyazeva et al., 1999 and unpublished), individual differences were noticed in the degree of interhemispheric synchronization. Their origin is at present unexplained. The two most likely possibilities are: differences in the morphology of callosal connections or individual variability in attention to the stimuli. Both possibilities could be tested in future experiments. The stimulus-induced differences emerge even more clearly from the population analysis shown in Fig. 9. The IG stimulus evokes synchronous activity over a cluster of electrodes overlying the occipital poles and spanning both hemispheres. On the other hand, the OG stimuli induce synchronized activity in two symmetrical clusters of electrodes overlying the parietotemporal regions of the two hemispheres. Statistical comparisons (permuted Hotelling’s T 2 test) returned highly significant differences (Bonferroni corrected for multiple comparisons; all electrodes included) for the OG and IG stimuli vs. background (P b 0.001) and for IG vs. OG (P b 0.005). At individual level, significant differences (P b 0.01) were obtained in all subjects in case the comparisons were restricted to the occipital electrodes with DS N 0.

Frequency bands arrangement Although the above analysis might be satisfactory within the more general frame of system theory, the prevailing trend in studies of oscillatory cortical activity is to analyze separate frequency bands. Hence, we broke down the analysis of the S estimator into the three traditional bands of EEG frequencies: the a (7–13 Hz), b (13–30 Hz), and c (30–70 Hz) frequency bands. This was also done to compare the results of the present method for detecting synchronization with those of coherence analysis (of the same EEG recordings) reported in (Knyazeva et al., 2002, 2003) and summarized in the next section. The results of such analysis at the population level are reported in Fig. 10. The bi-hemispheric, occipital synchronization induced specifically by the IG stimulus vs. background involves mainly the b band. The c band is also synchronized by the stimuli, but no differences are observed between the OG and the IG stimuli. Interestingly, the topography of synchronization in the b and c bands is also different. The selective activation in the b band is in agreement with the changes induced by the same stimuli in interhemispheric coherence (Knyazeva et al., 1999, 2002, 2003). However, surprisingly, the OG stimuli evoked two symmetrical clusters of synchronization in the parieto-temporal electrodes in the a band, and a weaker cluster of trans-hemispheric synchronization in

350

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

Fig. 8. Topographical maps of the variations of the mean synchronization (DS, see text) of EEG activity induced by the IG grating with respect to the OG grating. Other conventions as in Fig. 6.

the occipital lobe in the b band. They also evoked a cluster of transhemispheric occipito-parietal synchronization in the c band, similar to that induced by the IG stimulus.

Discussion Techniques for assessing synchronization of EEG activity, derived from system theory, have been recently proposed. Methodo-

Fig. 9. Average topographical maps, for the entire population of subjects, of the variations, with respect to the background, of the mean synchronization (DS, see text) of EEG activity induced either by the iso-oriented grating (IG vs. BGR) or by the orthogonally oriented grating (OG vs. BGR). Other conventions as in Fig. 6.

logically, these techniques can be categorized according to two main streams: those based on phase relationships (as reported for instance by Tass et al., 2003) and those based on state-space interdependences (e.g., Quiroga et al., 2002). Both have been employed on three fronts: firstly, at the brain modeling level, for verifying the hypothesis that synchronization exists between brain regions (Breakspear, 2002; Stam et al., 2003; Tononi et al., 1994); secondly, for investigating the functional role of synchronization within the neuronal populations involved in a task (Lachaux et al., 1999; Micheloyannis et al., 2003; Tass et al., 2003); and thirdly, in clinical environment, for assessing the effects of lesions or diseases in the brain (Quiroga et al., 2002; Mormann et al., 2000; Tass et al., 1998). Along the second methodological line, we have applied a method, based on the variation of the embedding state-space of a complex dynamical system, in order to investigate the functional role of a stimulus-induced EEG synchronization. Starting from the hypothesis that nonlinearities are an intrinsic feature in brain functioning, we introduced our measure for studying synchronization phenomena in conditions where linear methods (spectral coherence) have been already applied, which allows comparisons between the two methodologies. Our method is simpler than other nonlinear measures of synchronization (Boccaletti et al., 2002; Quiroga et al., 2002) and at least equally robust, making it suitable for the analysis of stimulus-driven EEG conditions,

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

351

Fig. 10. Average topographical maps, for the entire population of subjects, of the variations, with respect to the background, of the mean synchronization (S, see text) of EEG activity, within the a, b, and c frequency bands, induced either by the iso-oriented grating (IG vs. BGR) or by the orthogonally oriented grating (OG vs. BGR). Other conventions as in Fig. 6.

where other complex nonlinear approaches may turn out to be too fragile (Arnhold et al., 1999). Indeed, the S estimator used by this method allowed demonstrating stimulus-specific changes in the overall synchronization of brain activity, in the population analysis and partially at the level of the examined single subjects. We have applied S to simulated coupled oscillators, including situations caricaturing the hypothetical coupling generated between the hemispheres by IG and OG visual stimuli. In this test, S turned out to be at least as reliable and robust as the other estimators of synchronization. More importantly, S is naturally multivariate, which makes it particularly suitable for high-density surface sampling EEG analysis, currently replacing conventional EEG methods, both in research laboratories and in clinics. In the topographical part of the simulation, the coupling scheme caricaturing the OG condition returned two adjacent, but separate clusters of increased S. The scheme caricaturing the IG condition returned a wide continuous cluster extending across the hemispheres, comparable to the findings in the human EEG recordings. It is to be noticed though, that the displacement of the clusters driven by the different stimuli in the spatial and frequency domains was not grasped by the simulation. This has potentially interesting consequences for the interpretation of the EEG data, beyond the aims of the present study. The results obtained with the S estimator on EEG recordings can also be compared with those obtained with coherence analysis on the same material (Knyazeva et al., 2002, 2003). Because of its intrinsic multivariate nature addressing the whole skull topography, the S estimator turns out to be a complement to coherence analysis; the most commonly used linear method for the assessment of EEG synchronization. Indeed, the S estimator is based on the Principal Component Analysis (PCA), while coherence is based on the Discrete Fourier Transform (DFT).

However, it is proven that, considering a very large delay embedding state-space4, the PCA transformation tends to coincide with the DFT (Broomhead and King, 1986). Nevertheless, for small embedding dimension bthe PCA selects the relevant structures in [state] space rather than the relevant frequencies,Q (Kantz and Schreiber, 1997) and corresponds, therefore, to a synchronization analysis within the more general frame of system theory. Furthermore, when considering synchronization between more than a single pair of electrodes, the S estimator is more natural, being intrinsically multivariate, than coherence, and allows the analysis of synchronization phenomena spanning a wide territory. Finally, the S estimator addresses the investigations of broad-band synchronization phenomena. For the data considered in this study, the following differences between the results obtained with the S estimator and the EEG coherence analysis are noteworthy. The S estimator highlighted changes over a broader set of electrodes than the coherence analysis. The latter, applied to regions of interest, demonstrated changes mainly restricted to a single couple of electrodes in the occipital region, the electrodes 70 and 90 of Fig. 2, (Knyazeva et al., 2002, 2003). Significant coherence changes were restricted to a narrow band of frequencies between 20 and 30 Hz. Finally, the formation of lateralized clusters of synchronous activity in the temporal region in the a band was not seen with the coherence analysis, nor was the fact that synchronization in the c band did not differentiate between the two types of gratings. Theoretically,

4

For ergodic signals this approximately corresponds to very long epochs. Indeed, PCA is a diagonalization of a symmetric operator, and when the operator is also translationally-invariant (very long epochs), the eigenfunctions turn out to be complex exponentials, that is, the Fourier basis.

352

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

the results obtained with the S estimator might also be obtained with EEG coherence, though that would be computationally expensive. Indeed, because of the inherent bivariate nature of coherence, EEG coherence analysis is usually limited to couples of electrodes within a region of interest. In conclusion, both the EEG coherence study and the S estimator are useful tools for identifying topographical changes in the brain activity, which are the probable counterpart of functionally connected neuronal assemblies (for the definition of functional connectivity see Friston, 1994; Friston et al., 1993; Horwitz, 2003). The results obtained with the S estimator emphasize the complexity of the changes induced even by very simple stimuli, and appear to provide an easy to use and robust tool for the analysis of these changes. Although the S-based analysis might be satisfactory within the frame of system theory, from other viewpoints it is not. The prevailing concept of cortical functioning suggests that frequencyspecific oscillations are directly related to specific cortical functions, particularly in the visual domain (Eckhorn et al., 1988; Gray et al., 1989; Rodriguez et al., 1999). It would also be expected that the function- and frequency-specific oscillators are also localized at specific sites, whose identification is a challenging task for contemporary research. Therefore, in the second part of the study, we broke down the analysis of the S estimator to the a, b, and c frequency bands. This analysis demonstrated that, compared to the background, the identical gratings induced a coherent activation of the two hemispheres in the b band, focused on the occipital electrodes. Similar but weaker effects were obtained with the orthogonal stimuli. These results are in line with those obtained by coherence analysis. Unexpectedly, however, the orthogonal stimuli gave raise to two symmetrical clusters of activation in the a frequency range in the occipito-temporal region. The two clusters dissolved with the identical gratings. The c range was insensitive to stimulus diversity in the sense that both the orthogonal and collinear stimuli increased the values of the S estimator mostly in the parietal region. From the neurobiological point of view, the results summarized above raise two distinct sets of questions. The first one is related to the morpho-functional substrate of the observed phenomena. The response here is, to some extent, of an inferential kind. As one would have predicted by elementary notions of brain functional anatomy, the changes were concentrated in the occipital–temporal part of the brain. Indeed, the changes were induced by visual stimuli, which are known to powerfully activate the primary and secondary visual areas (Kiper et al., 1999; Knyazeva et al., 1999; Law et al., 1988). When the two hemispheres viewed identical gratings, a field of synchronized activity was formed, spanning across the hemispheres. Instead, when they viewed orthogonally oriented and moving gratings, two separate fields of synchronized activity were formed in the two hemispheres. Presumably, the recorded activities inform on events occurring in relative close proximity to the electrodes involved, that is, in the occipital temporal region. More precision on the topographical location of the recorded signal can be derived only from techniques of dipole localization (Gevins et al., 1995; Sidman et al., 1991), or possibly by correlation with the fMRI activation (Knyazeva et al., 2003). Concerning the substrate of the interhemispheric synchronization, the similarity of some of our findings with those obtained by coherence responses to the same stimuli in animals and in man (Kiper et al., 1999; Knyazeva et al., 1999), suggests that, as in the latter, the cortico-

cortical connections mediated through the corpus callosum, might be crucially involved. The second question is of a more general, theoretical, and methodological nature. The observed phenomena could be explained as the result of the interaction between different classes of cortical oscillators, each class being characterized by a different oscillatory frequency (Breakspear, 2002; Corchs and Deco, 2001; Frank et al., 2000). Different stimuli could differently bind/unbind the oscillators. Alternatively, the stimuli themselves could be involved in the generation of oscillations at different frequencies as well as in their binding. Possibly, some of the uncertainties discussed here could be resolved by the systematic analysis of periodic visual stimuli of different frequencies, and/or of non periodic visual stimuli. Irrespectively of the answers to the questions above, another consideration seems to emerge from the results we described. The response of the visual areas to the extremely simple stimuli used is neither simple nor localized in either the frequency or in the spatial domain. Hence, as already suspected (Basar et al., 2001; Bhattacharya and Petsche, 2002; Bressler et al., 1993; Chen et al., 2003), it is unlikely that the response might be satisfactorily characterized by assuming the existence of simple relations between EEG frequencies and functions. Instead, the complexity of cortical dynamics might be better characterized by approaches considering that brain activity continuously varies with time within a complex state-space (Bruns et al., 2000; Eckhorn et al., 2001). Such a method might be derived from application of chaotic models as several research lines in the literature seem to suggest (Dafilis et al., 2001; Freeman, 2003; Kaneko and Tsuda, 2003; Kay, 2003; Wright and Liley, 1996).

Conclusions We introduced S, a new estimator that indirectly quantifies the amount of synchronization within a data set measuring the contraction of the embedding dimension of the state-space, which is a theoretical consequence of synchronization. The new estimator is grounded on the Principal Components Analysis embedding technique. Applications to simulated coupled oscillators and to EEG data indicated that S provides a robust new tool for the analysis of synchronous neuronal assemblies in the cerebral cortex. In the current study, the results obtained with this tool were comparable to, albeit richer than those obtained with the more traditional method of EEG coherence analysis. The S estimator revealed that very simple visual stimuli may cause complex changes in the dynamics of cortical activity and emphasized the need of analytical tools to study brain functioning from a holistic viewpoint, albeit reducible to the dynamics of system components. Finally, because of its intrinsic multivariate nature, this new estimator is particularly suitable for multichannel EEG analysis, and is a good candidate for clinical applications, such as online analysis and monitoring of the state of functional connectivity.

Acknowledgments We are grateful to Dr. Ute Feldmann for the exciting discussions and for having inspired some of the ideas, and an anonymous reviewer for his helpful comments on earlier drafts.

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

This work was supported by the European project APEREST: IST2001-34893 and OFES-01.0456 and by the Swiss National Foundation grant 31-63894.00.

References Anokhin, A.P., Lutzenberger, W., Birbaumer, N., 1999. Spatiotemporal organization of brain dynamics and intelligence: an EEG study in adolescents. Int. J. Psychophysiol. 33, 259 – 273. Aoki, F., Fetz, E.E., Shupe, L., Lettich, E., Ojemann, G.A., 1999. Increased gamma-range activity in human sensorimotor cortex during performance of visuomotor tasks. Clin. Neurophysiol. 110, 524 – 537. Arnhold, J., Grassberger, P., Lehnertz, K., Elger, C.E., 1999. A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D 134, 419 – 430. Basar, E., Basar-Eroglu, C., Karakas, S., Schurmann, M., 2001. Gamma, alpha, delta and theta oscillations govern cognitive processes. Int. J. Psychophysiol. 39, 241 – 248. Bhattacharya, J., Petsche, H., 2002. Shadows of artistry: cortical synchrony during perception and imagery of visual art. Cognit. Brain Res. 13, 179 – 186. Boccaletti, S., Kurths, J., Osipov, G., Valladares, D.L., Zhou, C.S., 2002. The synchronization of chaotic systems. Phys. Rep. 366, 1 – 101. Breakspear, M., 2002. Nonlinear phase desynchronization in human electroencephalographic data. Hum. Brain Mapp. 15, 175 – 198. Bressler, S.L., Coppola, R., Nakamura, R., 1993. Episodic multiregional cortical coherence at multiple frequencies during visual task performance. Nature 366, 153 – 156. Brown, R., Kocarev, L., 2000. A unifying definition of synchronization for dynamical systems. Chaos 10, 344 – 349. Bruns, A., Eckhorn, R., Jokeit, H., Ebner, A., 2000. Amplitude envelope correlation detects coupling among incoherent brain signals. NeuroReport 11, 1509 – 1514. Broomhead, D.S., King, G.P., 1986. Extracting qualitative dynamics from experimental data. Physica D 20, 217 – 236. Burgess, A.P., Gruzelier, J.H., 1997. Short duration synchronization of human theta rhythm during recognition memory. NeuroReport 8, 1039 – 1042. Corchs, S., Deco, G., 2001. A neurodynamical model for selective visual attention using oscillators. Neural Netw. 14, 981 – 990. Chen, Y., Ding, M., Kelso, J.K., 2003. Task-related power and coherence changes in neuromagnetic activity during visuomotor coordination. Exp. Brain Res. 148, 105 – 116. Dafilis, M.P., Liley, D.T.J., Cadusch, P.J., 2001. Robust chaos in a model of the electroencephalogram: implications for brain dynamics. Chaos 11, 474 – 478. De Feo, O., Maggio, G.M., Kennedy, M.P., 2000. The Colpitts oscillator: families of periodic solutions and their bifurcations. Int. J. Bifurc. Chaos 10, 935 – 958. Eckhorn, R., Bauer, R., Jordan, W., Brosch, M., Kruse, W., Munk, M., Reitboeck, H.J., 1988. Coherent oscillations: a mechanism of feature linking in the visual cortex? Multiple electrode and correlation analyses in the cat. Biol. Cybern. 60, 121 – 130. Eckhorn, R., Bruns, A., Saam, M., Gail, A., Gabriel, A., Brinksmejer, H.J., 2001. Flexible cortical gamma-band correlations suggest neural principles of visual processing. Vis. Cogn. 8, 519 – 530. Engel, A.K., Konig, P., Kreiter, A.K., Singer, W., 1991. Interhemispheric synchronization of oscillatory neuronal responses in cat visual cortex. Nature 252, 1177 – 1179. Fell, J., Fernandez, G., Elger, C.E., 2003. More than synchrony: EEG chaoticity may be necessary for conscious brain functioning. Med. Hypotheses 61, 158 – 160. Ford, M.R., Goethe, J.W., Dekker, D.K., 1986. EEG coherence and power changes during a continuous movement task. Int. J. Psychophysiol. 4, 99 – 110.

353

Frank, T.D., Daffertshofer, A., Peper, C.E., Beek, P.J., Haken, H., 2000. Towards a comprehensive theory of brain activity: coupled oscillator systems under external forces. Physica D 144, 62 – 86. Freeman, W.J., 2003. Evidence from human scalp electroencephalograms of global chaotic itinerancy. Chaos 13, 1067 – 1077. Friston, K.J., 1994. Functional and effective connectivity in neuroimaging: a synthesis. Hum. Brain Mapp. 2, 56 – 78. Friston, K.J., Frith, C.D., Liddle, P.F., Frackowiak, R.S.J., 1993. Functional connectivity: the principal-component analysis of large (PET) data sets. J. Cereb. Blood Flow Metab. 13, 5 – 14. Gevins, A., Leong, H., Smith, M.E., Le, J., Du, R., 1995. Mapping cognitive brain function with modern high-resolution electroencephalography. Trends Neurosci. 18, 429 – 436. Goldman, R.I., Stern, J.M., Engel, J.J., Cohen, M.S., 2000. Acquiring simultaneous EEG and functional MRI. Clin. Neurophysiol. 111, 1974 – 1980. Gray, C.M., 1999. The temporal correlation hypothesis of visual feature integration: still alive and well. Neuron 24, 31 – 47. Gray, C.M., Kfnig, P., Engel, A.K., Singer, W., 1989. Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature 338, 334 – 337. Higgins, J.J., 2004. Introduction to Modern Nonparametric Statistics. Brooks/Cole-Thomson Learning, USA. Horwitz, B., 2003. The elusive concept of brain connectivity. NeuroImage 19, 466 – 470. Inouye, T., Shinosaki, K., Sakamoto, H., Toi, S., Ukai, S., Iyama, A., Katzuda, Y., Hirano, M., 1991. Quantification of EEG irregularity by use of the entropy of the power spectrum. Electroencephalogr. Clin. Neurophysiol. 79, 204 – 210. Johnson, R., Wichern, D., 2002. Applied Multivariate Statistical Analysis. Prentice Hall, Englewood Cliffs, NJ, USA. Joliffe, I., 2002. Principal Component Analysis, 2nd ed. Springer-Verlag, New York, NY, USA. Kaneko, K., Tsuda, I., 2003. Chaotic itinerancy. Chaos 13, 926 – 936. Kantz, H., Schreiber, T., 1997. Nonlinear Time Series Analysis. Cambridge Univ. Press, Cambridge, UK. Kay, L.M., 2003. A challenge to chaotic itinerancy from brain dynamics. Chaos 13, 1057 – 1066. Keil, A., Muller, M.M., Ray, W.J., Gruber, T., Elbert, T., 1999. Human gamma band activity and perception of a gestalt. J. Neurosci. 19, 7152 – 7161. Kiper, D., Knyazeva, M., Tettoni, L., Innocenti, G., 1999. Visual stimulusdependent changes in interhemispheric EEG coherence in ferrets. J. Neurophysiol. 82, 3082 – 3094. Klimesch, W., 1999. EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29, 169 – 195. Klimesch, W., Schimke, H., Schwaiger, J., 1994. Episodic and semantic memory: an analysis in the EEG theta and alpha band. Electroencephalogr. Clin. Neurophysiol. 91, 428 – 441. Knyazeva, M., Kiper, D., Vildavski, V., Despland, P., Maeder-Ingvar, M., Innocenti, G., 1999. Visual stimulus-dependent changes in interhemispheric EEG coherence in humans. J. Neurophysiol. 82, 3095 – 3107. Knyazeva, M.G., Maeder, P., Fornari, E., Meuli, R., Innocenti, G.M., 2002. Interhemispheric synchronization under visual stimulation: the evidence from dense array EEG. Clin. Neurophysiol. 113 (Suppl. 1), S68 (abstract). Knyazeva, M.G., Fornari, E., Meuli, R., Innocenti, G.M., Maeder, P., 2003. Interhemispheric EEG synchronization correlates with fMRI BOLD response. International Conference on Human Brain Mapping, June 18– 22, 2003, New York, USA, NeuroImage, vol. 19, p. N2. Available on CD-ROM in, (abstr. 1478). Koeda, T., Knyazeva, M., Njiokiktjien, C., Jonkman, E.J., Sonneville, L.D., Vildavsky, V., 1995. The EEG in acallosal children. coherence values in the resting state: left hemisphere compensatory mechanism? Electroencephalogr. Clin. Neurophysiol. 95, 397 – 407. Krishnamoorthy, K., Jianqi, Y., 2004. Modified Nel and Van der Merwe test

354

C. Carmeli et al. / NeuroImage 25 (2005) 339–354

for the multivariate Behrens-Fisher problem. Stat. Probab. Lett. 66, 161 – 169. Kuks, J.B., Vos, J.E., O’Brien, M.J., 1987. Coherence patterns of the infant sleep EEG in absence of the corpus callosum. Electroencephalogr. Clin. Neurophysiol. 66, 8 – 14. Lachaux, J.P., Rodriguez, E., Martinerie, J., Varela, F.J., 1999. Measuring phase synchrony in brain signals. Hum. Brain Mapp. 8, 194 – 208. Law, M.I., Zahs, K.R., Stryker, M.P., 1988. Organization of primary visual cortex (Area 17) in the ferret. J. Comp. Neurol. 278, 157 – 180. Lemieux, L., Krakow, K., Fish, D.R., 2001. Comparison of spike-triggered functional MRI BOLD activation and EEG dipole model localization. NeuroImage 14, 1097 – 1104. Lopes Da Silva, F.H., Vos, J.E., Mooibroek, J., Van Rotterdam, A., 1980. Partial coherence analysis of thalamic and cortical alpha rhythms in dog—A contribution towards a general model of the cortical organization of rhythmic activity. In: Pfurtscheller, G., Buser, P., Lopes da Silva, F.H., Petsche, H. (Eds.), Rhythmic EEG Activities and Cortical Functioning. Elsevier, Amsterdam, The Netherlands. Lorenz, E.N., 1963. Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130 – 141. Micheloyannis, S., Vourkas, M., Bizas, M., Simos, P., Stam, C.J., 2003. Changes in linear and nonlinear EEG measures as a function of task complexity: evidence for local and distant signal synchronization. Brain Topogr. 15, 238 – 247. Miltner, W.H.R., Braun, C., Arnold, M., Witte, M., Taub, E., 1999. Coherence of gamma-band EEG activity as a basis for associative learning. Nature 397, 434 – 436. Mima, T., Oluwatimilehin, O., Hiraoka, T., Hallet, M., 2001. Transient interhemispheric neuronal synchrony correlates with object recognition. J. Neurosci. 21, 3942 – 3948. Montplaisir, J., Nielsen, T., Cote, J., Boinvin, D., Rouleau, I., Lapierre, G., 1990. Interhemispheric EEG coherence before and after partial callosotomy. Clin. Electroencephalogr. 21, 42 – 47. Mormann, F., Lehnertz, K., David, P., Elger, C.E., 2000. Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients. Physica D 144, 358 – 369. Mulert, C., Jager, L., Pogarell, O., Bussfeld, P., Schmitt, R., Juckel, G., Hegerl, U., 2002. Simultaneous ERP and event-related fMRI: focus on the time course of brain activity in target detection. Methods Find. Exp. Clin. Pharmacol. 24, 17 – 20. Muller, M.M., Bosch, J., Elbert, T., Kreiter, A., Sosa, M.V., Sosa, P.V., Rockstroh, B., 1996. Visually induced gamma-band responses in human electroencephalographic activity—A link to animal studies. Exp. Brain Res. 112, 96 – 102. Pinkofsky, H.B., Struve, F.A., Meyer, M.A., Patrick, G., Reeves, R.R., 1997. Decreased multi-band posterior interhemispheric coherence with a lipoma on the corpus callosum: a case report of a possible association. Clin. Electroencephalogr. 28, 155 – 159. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in C: the Art of Scientific Computing, 2nd ed. Cambridge Univ. Press, Cambridge, NY, USA. Quarteroni, A., Sacco, R., Saleri, F., 2000. Numerical Mathematics. Springer-Verlag, Berlin, Germany.

Quiroga, R.Q., Arnhold, J., Lehnertz, K., Grassberger, P., 2000. KullbackLeibler and renormalized entropies: applications to electroencephalograms of epilepsy patients. Phys. Rev., E 62, 8380 – 8386. Quiroga, R.Q., Kraskov, A., Kreuz, T., Grassberger, P., 2002. Performance of different synchronization measures in real data: a case study on electroencephalographic signals. Phys. Rev., E 65 (Art. 041903). Rodriguez, E., George, N., Lachaux, J., Martinerie, J., Renault, B., Varela, F., 1999. Perception’s shadow: long-distance synchronization of human brain activity. Nature 397, 430 – 433. Rfssler, O.E., 1976. Chaotic behavior in simple reaction systems. Z. Naturforsch. 31A, 259 – 264. Sabesan, S., Narayanan, K., Prasad, A., Spanias, A., Sackellares, J., Iasemidis, L., 2003. Predictability of epileptic seizures: a comparative study using Lyapunov exponent and entropy based measures. Biomed. Sci. Instrum. 39, 129 – 135. Schiff, S.J., So, P., Chang, T., Burke, R.E., Sauer, T., 1996. Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble. Phys. Rev., E 54, 6708 – 6724. Sidman, R.D., Major, D.J., Ford, M.R., Ramsey, G.G., Schlichting, C., 1991. Age-related features of the resting pattern-reversal visual evokedresponse using the dipole localization method and cortical imaging technique. J. Neurosci. Methods 37, 27 – 36. Srinivasan, R., 1999. Methods to improve the spatial resolution of EEG. Int. J. Bioelectromagn. 1, 102 – 111. Srinivasan, R., 2003. High-resolution EEG: theory and practice. In: Handy, T.C. (Ed.), Event Related Potentials: a Methods Handbook. MIT Press, Cambridge, MA. Stam, C.J., Breakspear, M., van Walsum, A.M.V., van Dijk, B.W., 2003. Nonlinear synchronization in EEG and whole-head MEG recordings of healthy subjects. Hum. Brain Mapp. 19, 63 – 78. Tallon, C., Bertrand, O., Bouchet, P., Pernier, J., 1995. Gamma-range activity evoked by coherent visual stimuli in humans. Eur. J. Neurosci. 7, 1285 – 1291. Tallon-Baudry, C., Bertrand, O., Fischer, C., 2001. Oscillatory synchrony between human extrastriate areas during visual short-term memory maintenance. J. Neurosci. 21 (1–5), RC177. Tass, P., Rosenblum, M.G., Weule, J., Kurths, J., Pikovsky, A., Volkmann, J., Schnitzler, A., Freund, H., 1998. Detection of n:m phase locking from noisy data: application to magnetoencephalography. Phys. Rev. Lett. 81, 3291 – 3294. Tass, P.A., Fieseler, T., Dammers, J., Dolan, K., Morosan, P., Majtanik, M., Boers, F., Muren, A., Zilles, K., Fink, G.R., 2003. Synchronization tomography: a method for three-dimensional localization of phase synchronized neuronal populations in the human brain using magnetoencephalography. Phys. Rev. Lett. 90 (Art. 088101). Tononi, G., Sporns, O., Edelman, G.M., 1994. A measure for brain complexity: relating functional segregation and integration in the nervous system. Proc. Natl. Acad. Sci. U. S. A. 91, 5033 – 5037. Wackermann, J., 1996. Beyond mapping: estimating complexity of multichannel EEG recordings. Acta Neurobiol. Exp. 56, 197 – 208. Wright, J.J., Liley, D.T.J., 1996. Dynamics of the brain at global and microscopic scales: Neural networks and the EEG. Behav. Brain Sci. 19, 285 – 320.