Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs.

Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs. J The analytic signal given by the Hilbert transform applied to an electroencephalographic (EEG) trace is a vector of instan- taneous amplitude and phase at the temporal resolution of the digitizing interval (here 2 ms). The transform was applied after band-pass ﬁltering for extracting the gamma band (20–80 Hz in rabbits) to time series from up to 64 EEG channels recorded simultaneously from high-density arrays giving spatial “windows” of 4 (cid:1) 4 to 6 (cid:1) 6 mm onto the visual, auditory, or somatosensory cortical surface. The time series of the analytic phase revealed phase locking for brief time segments in spatial patterns of nonzero phase values from multiple EEG that was punctuated by episodic phase decoherence. The derivative of the analytic phase revealed spikes occurring not quite simul- taneously (within (cid:2) 4 ms) across arrays aperiodically at mean rates in and below the theta range (3–7 Hz). Two measures of global syn- chronization over a group of channels were derived from analytic phase differences between pairs of channels on the same area of cortex. One was a synchronization index expressing phase locking. The other was a decoherence index estimating the variance in phase among multiple channels. Spectral analyses of the indices indicated that decoherence events recurred aperiodically at rates in and below the theta range of the EEGs. The results provide support for the hypothesis that neurons in mesoscopic neighborhoods in sensory cortices self-organize their activity by synaptic interactions into wave packets that have spatial patterns of amplitude (AM) and phase (PM) modulation of their spatially coherent carrier waves in the gamma range and that form and dissolve aperiodically at rates in and below the theta range. Each AM pattern is formed by a nonlinear state transition in the cortical dynamics, as shown by spikes in the deriv- ative. Phase locking within each PM pattern is not at zero phase lag but over a ﬁxed distribution of phase values that is consistent with the radially symmetric phase gradients already reported called “phase cones” detected by Fourier-based methods. The insight is suggested that sensory cortices are bistable comparably to cardiac dynamics, with a diastolic state that accepts sensory input and an abrupt transi- tion to a systolic state that transmits perceptual output. Further support for this inference will require improvements in methods for temporal resolution of the times of onset of spatial patterns of phase modulation.


I N T R O D U C T I O N
Episodic synchronization among selected neuronal pulse trains in local areas of cortex has suggested phase locking of neural activity in response to sensory stimuli (e.g., Singer and Gray 1995). In support, EEGs from primary sensory areas with high-density electrode arrays in rabbits and cats have revealed repeated episodes of spatially coherent gamma oscillations (Barrie et al. 1996;Freeman and Viana Di Prisco 1986) called "wave packets" (Freeman 1975). These packets were observed through 4 ϫ 4 to 6 ϫ 6 mm "windows" provided by epidural 8 ϫ 8 electrode arrays with interelectrode spacing of 0.5-0.8 mm fixed over each of the sensory areas prior to classical or operant conditioning. Up to 64 simultaneously recorded samples of the aperiodic carrier wave of the wave packets showed spatial patterns of amplitude modulation (AM) that served to classify brief electroencephalographic (EEG) epochs with respect to olfactory, visual, auditory, or somatosensory conditioned stimuli. Wave packets in rabbits were approximately 80 -100 ms in duration (Freeman and Baird 1987;Freeman and Barrie 2000), and they recurred irregularly at mean rates in the theta range (2-7 Hz). Wave packets were accompanied by spatial patterns of phase modulation (PM) in the form of a cone, for which the diameter ranged up to 3 cm or more, and for which the location and sign of the apex varied randomly with each successive wave packet, without relation to conditioned stimuli or anatomical landmarks.
The inference was made from the breadth of nonzero synchrony revealed by the phase values that the wave packets formed by cooperative synaptic interactions within areas of cortex on a mesoscopic scale greater than the microscopic space-time scale of the component neurons, but smaller than the macroscopic scale of domains that have been shown by brain imaging based on metabolic and hemodynamic processes (Freeman 2000a). These mesoscopic interactions shaped the spatial AM patterns in accordance with widespread changes in intracortical synaptic strengths from cumulative effects of reinforcement learning.
Each successive wave packet with its spatial AM pattern was formed by a nonlinear state transition having a site of nucleation demarcated by the apex of a conical phase gradient and a radial velocity determined by the conduction velocities of intracortical axons running parallel to the pial surface. The cumulative phase difference with distance from the apex was conceived to provide a soft boundary condition limiting the size of the mesoscopic cooperative domains to the half-power radius at Ϯ45°, thus giving an estimated state transition time on the order of 3-7 ms and a modal diameter on the order of 5-15 mm, considerably greater than the dimensions of the dendritic arbors of most neurons, the estimated sizes of cortical columns and hypercolumns, and the size of the array window, but less than the dimensions of the lobes of the forebrain and, in small mammals, the cerebral hemispheres (Freeman 2000d;Freeman and Barrie 2000).
The significance of these findings and inferences lies in the possibility that the episodes of spatial synchrony manifest a mesoscopic integrative process by which information that is injected into cortical activity by sensory receptors through thalamic relays is integrated into perceptual forms shaped by life-long prior learning. If so, the process offers an alternative to a widely debated linear form of the "binding" hypothesis (Singer and Gray 1995;Tovée and Rolls 1992;von der Malsburg 1983), which holds that the pulse trains of a sparse network comprising a small number of neurons embedded in cortical neuropil are locked in oscillation at a single frequency with zero time lag synchrony (e.g., Traub et al. 1996). The alternative is that the binding is at the mesoscopic level and consists of a small covariant fraction of the total variance of all the neurons within the diameter of a wave packet having a common, broad-spectrum oscillation (Freeman 2000a). The covariance imposed by cooperative dynamics gives a shared wave form in the gamma range despite its lack of confinement to one frequency. Phase values can be defined at the peak or center frequency in the spectrum of the oscillations over the duration of each wave packet with respect to the phase of the spatial ensemble average. Clearly the cooperative activity depends on communication by action potentials among cortical neurons with requisite propagation delays, which impose a distribution of phase values over its spatial extent of cooperation (Freeman 2000c), leading to the appearance of phase cones. Phase differences among components in networks having periodic oscillations can be steady for those with a fixed frequency, but, in chaotic systems and especially those having noise, the phase differences between components are not steady. However, if there is significant strength in the cooperativity, the components may show plateaus of steady phase differences surrounded by episodes of increasing differences known as "phase slip" (Pikovsky et al. 1997) that may resemble a random walk process and that reflect changing levels of coupling or noise or both.
Pursuit of this alternative hypothesis (Freeman 2000d) required a new approach with better temporal resolution of the widespread changes in phase at the times of state transitions when new phase cones formed. Finer grain of measurement was given by the Hilbert transform, by which each EEG trace is expressed in two time series: the instantaneous analytic amplitude and analytic phase at each digitizing step. The results reported here complement the measurements of phase using the Fourier transform on the same data, with sufficient resolution to demonstrate that abrupt phase jumps may occur in each trace in less than the 2-ms digitizing step, and not simultaneously in multiple traces but over a range of time intervals across electrode arrays, as predicted by prior estimates from measurements of phase cones of the time required for state transitions to spread spatially from sites of nucleation (Freeman 2000d;Freeman and Barrie 2000).

M E T H O D S
The methods have been described in prior publications by which the EEG data used here were acquired from rabbits (Barrie et al. 1996;Freeman and Barrie 2000) and cats (Gaál and Freeman 1997), including the techniques of animal surgery, veterinary care, design and assembly of electrode arrays, training to olfactory, visual, auditory, and somatosensory conditioned stimuli (CS) with reinforcement (CSϩ) or without (CSϪ), recording of up to 64 EEG traces simultaneously in blocks of 40 trials each lasting 6 s, and data preprocessing and storage. The conditioned response (CR) in rabbits was the autoshaped response of sniffing as a part of the orienting reflex (Freeman et al. 1983) in a classical aversive paradigm with weak electric shock to the cheek for reinforcement. The CR in cats was a bar press for food in an operant appetitive paradigm.
Likewise the basic theory and digital methods for computing the Hilbert transform have been well described (Barlow 1993;Hahn 1996;Stearns and David 1996). Essentially a fluctuating time series is decomposed into two time series, one showing the analytic amplitude, v(t), and the other the analytic phase, (t), both as functions of time.
Branching of the arc-tangent occurs whenever the signal crosses zero, giving a sawtooth function, so the phase must be "unwrapped" to form a continuously increasing function of time. Matlab (Stearns and David 1996) provided an unwrap function that was not suitable for these data, because it sensed abrupt changes of the sawtooth series larger than a criterion and added a constant amount to correct. The analysis of phase cones showed that the step size of the apparent changes with state transitions might in theory have any value from Ϫ to ϩ without a threshold, although the observed steps would be smoothed by the effects of local volume conduction. The sawtooth changes at branch points requiring unwrapping also varied over a wide range. Setting a criterion to find them too low introduced false detection. However, the sawtooth changes only occurred near zero crossings of the filtered EEG, but the physiological jumps were not restricted in that way, so an unwrap function was devised that used the zero crossings of the signal to predict the branch points and add a correction according to the size of the sawtooth change.
Extracting the analytic phase from systems having broad-spectrum chaotic dynamics (Lai et al. 1999) and especially having noise (García-Ojalvo and Sancho 1999;Rosenblum et al. 1998;Tass et al. 1998) is attended by substantial difficulties owing to the high variance and unpredictability of the estimates of instantaneous phase under those conditions. Prior experimental studies (Barrie et al. 1996) and theoretical modeling (Freeman 2000a;Kozma and Freeman 2001) of the EEG data used in the present studies showed that the temporal spectra conform to "1/f ␣ " distribution with ␣ ϭ Ϫ1.96 Ϯ 0.31 (mean Ϯ SE), that spatial patterns of amplitude and phase relating to behavior were found in the gamma range (20 -80 Hz in rabbits; 35-60 Hz in cats), and that optimal classifications of the AM patterns were obtained after band-pass filtering in these ranges. The 6-s trial segments were filtered for gamma activity in these bands using a Hanning window. The advantage of filtering before using the Hilbert transform was that large low-frequency baseline shifts and small high-frequency fluctuations were less likely to obscure the locations of gamma-band zero-crossings. Phase slip was reduced but not removed, particularly in times with low-amplitude and relatively higher noise. After filtering, the analytic signals were calculated in 3-s time segments starting 1 s prior to the stimulus, using the Matlab (version 5) Hilbert transform. Unlike the outcome of the Fourier transform in which frequency was determined first and then phase, the Hilbert transform first gave sequential values of phase in which frequency was implicit. Average frequency in radians/s (2f in Hz) of single-channel EEGs over time segments of 3 s starting 1 s before CS onset was estimated by fitting a line to the phase values, (t) in radians, where ϭ 1, 64 designated the channel or trace (electrode) number. Instantaneous frequency was numerically approximated in radians/second from phase differences between each digitizing step and the one two bins (4 ms) earlier. The sequence corresponded to the phase derivative, d(t)/dt in radians/second, as a function of time for each channel.
Phase differences between pairs of channels in a rectangular grid were calculated by subtracting their time series over 3 s in each trial. Synchronization between a pair of channels was inferred from a statistical tendency to maintain a nearly constant phase difference over a given time period even as the analytic phase changed markedly with time. Global synchronization was defined as a shared tendency to constancy among all pairs of channels in a grid fixed on one cortical area. This definition suggested a straightforward measure of synchronization using the SD, SD x , of the phase differences, ⌬ x1,x2 (t), between two channels, x1 and x2, after the local means were taken into consideration. This was done by subtracting the mean for each channel, x, from the values of ,x (t) in an 80 ms-window (long enough to encompass several cycles of gamma activity) The temporal SD from that mean, SD t (t), was calculated for the entire set of the normalized values in the window starting at t for each trace and channel x The spatial SD was computed for the phase differences at each time step in the window

͑3)
The average was calculated of the 40 values over the time window The time series, SD x (t) and SD t (t), were given by repeating the calculations while stepping the window at 10-ms intervals. They were generalized to a global measure for a group of channels by averaging the variances across all of the channel pairs and taking the square root. The value of this measure lay in its independence from the mean of the phase differences that were maintained between all channel pairs, which varied across the groups (8 ϫ 8 electrodes fixed in a 6 ϫ 6 mm window on each sensory cortex in rabbits and up to 16 electrodes in sparser 4 ϫ 4 mm windows on 4 cortices in 4 cats and 3 rabbits).
Since SD x (t) varied inversely with synchronization, it was used as a decoherence index. The values of ⌬ , (t), SD t (t), or SD x (t) were pooled in histograms to determine their distributions, as the basis for devising levels of significance for indices of the phase variation coarse-grained by the 80-ms window. An alternative and superior pairwise index of synchronization was adapted from a method developed by Tass et al. (1998), who designed their index to evaluate coupling between pairs of nonlinear systems for magnetoencephalographic and muscle signals observed in patients with Parkinson's disease. They defined statistical phase locking as a peak in the distribution of the phase differences between pairs of traces within a sliding window, based on normalized Shannon entropy where p j was the relative frequency of finding the phase mod 2 within the jth bin, and e varied between zero and ê ϭ ln N, the number of bins (e.g., 100 bins of 0.06 rad between Ϫ and ϩ radians). This synchronization index as a function of time was normalized so that q(t) was zero for a uniform distribution, and one for a delta distribution of phase values. They subtracted the 95.4th percentile of q calculated for a simulated data set composed of white noise, filtered in the same way as the data set, to establish a statistical significance level. This step was omitted here in favor of alternative controls. One control analyzed white noise treated in the way the EEG signals were, as used by Tass et al. (1998). Significance was inferred by comparing the values derived from the experimental data with results from white noise (time series of normally distributed random numbers with zero mean and unit SD, filtered, truncated to 12 bits, and otherwise treated in the same way as the EEG data). A second control was the use of "shuffled" time series, in which the start of each window was randomly reset independently on each channel, and the beginning segment was transposed and attached to the end. Thus each channel's phase time series of n points in a time window, (t), was replaced by (1, 2, . . . , k Ϫ 1)], and k was a point chosen from a uniformly distributed random variable with values from 1 to n. This shuffling in time guaranteed the independence of the channels.
The synchronization index, q(t), was generalized to multiple channels by calculating the distribution of phase differences over all pairs of channels, after subtracting the means for each pair within the 80-ms sliding window as in Eq. 3. This procedure imposed near-zero mean on all pairs, to document the phase locking of multiple traces that were fluctuating together in the gamma range but with phase differences that were transiently fixed at different values for every pair of traces, as has been found in the form of phase cones in all primary sensory receiving areas. plitude, v(t), in arbitrary units (top dashed curve), and the analytic phase, (t), in radians (sawtooth curve). Figure 2 shows the monotonically increasing values of (t) in radians after "unwrapping" for a selection of EEG traces from the visual and entorhinal cortices. The slopes reflect the average frequencies (near 250 rad/s ϭ 40 Hz) over the 3-s display in which a stimulus arrived at 0 ms. Cumulative small differences in phase between channels led to differing slopes and mean frequencies. Figure 3 (top frame) shows the successive values of the spatial differences in phase, ⌬ x1,x2 (t), between representative pairs of channels for every time point t ϭ Ϫ1,000 to ϩ2,000 ms from visual cortical EEGs. The "staircase" appearance indicated that phase differences were relatively constant over time (phase locking) between brief episodes of rapid change in either direction and not rising monotonically on any channel.
The fine structure of the decoherence episodes was examined by taking the numerical derivatives of the analytic phase values, d/dt. These derivatives highlighted the abruptness of the phase changes that occurred on multiple channels of intracortical EEGs (Fig. 3, bottom frame). The brief durations of the spikes in phase derivative values on each trace, and the variation in magnitude of the spikes indicated that the time taken by the jump at each channel might be below the 2-ms resolution of digitizing, although the appearance in time might be spread out by volume conduction. The spread of Ϯ2-3 bins (Ϯ4 -6 ms) in the times of detection was consistent with the known phase gradients revealed by phase cones in these data. The variation in magnitude was consistent with the replacement of one phase cone by another having a different location and sign of apex at each decoherence. The smoothness of the temporal change in phase difference between each pair was consistent with the effects of volume conduction of the local field potentials.
Global decoherence in the phase derivative time series was identified by two criteria. First, the height on a channel had to exceed 1 SD of the d/dt values (in practice 0.3 rad/4 ms) measured relative to the local mean surrounding the peak on that channel. Second, based on a histogram showing how many spikes by the first criterion occurred on multiple channels as a function of the number of channels, the criterion was adopted that spikes had to be detected on at least 60% of the channels within Ϯ4 ms. The larger tick marks on the abscissa of the bottom frame in Fig. 3 show the times at which widespread decoherence was detected with these two criteria.
Display of the high temporal resolution showed that the step changes were not instantaneous at each site or co-temporal across sites but spread over time intervals between 4 and 10 ms, and with the size of differences ranging in both directions from zero between different pairs up to 12 rad or more. An example in Fig. 4 was taken from the jump in Fig. 3 between 1,280 and 1,300 ms. The analysis showed that the jumps were not owing to failures in the methods used for unwrapping, because the phase changes took place over several milliseconds. FIG. 2. Examples are shown of the analytic phase, (t), for a representative set of 8 traces from the visual and entorhinal cortices of a cat. After unwrapping, (t) increased monotonically with time at a rate close to 250 rad/s (40 Hz). The differences in slope of traces recorded simultaneously from the same cortex often exceeded those between cortices.
FIG. 3. Top frame: differences in phase for a given pair of channels, ⌬ x1,x2 (t) ϭ x2 (t) Ϫ x1 (t), are plotted for several channel pairs, for one trial (cat data). These showed emergent phase differences between channels having the form of plateaus punctuated by abrupt changes nearly simultaneously on multiple channels. The direction of the jump varied at random irrespective of prior jumps, so the differences did not increase monotonically. This pattern was predicted from prior measurement of phase cones, in which the location and sign of the apices varied randomly (Freeman and Baird 1987;Freeman and Barrie 2000) and so also the direction of the phase gradient between each pair of electrodes. Bottom frame: the derivatives of the analytic phase, d/dt, were calculated for the above traces and plotted as time series revealing irregularly spaced spikes on multiple channels, which varied markedly in direction and magnitude. The variations could be attributed to the facts that the duration of the jumps in phase was close to the digitizing interval, and the direction of the change in spatial phase gradient had been found to change unpredictably between pairs of electrodes. The tick marks on the abscissa show the times of occurrence of global transitions in the sets of up to 64 traces in accord with criteria described in the text. A summary histogram of the intervals between them is shown in Fig. 8, and a power density spectrum is shown in Fig. 9.
Histograms of the distribution of ⌬ x1,x2 (t) are shown in Fig.  5, which provided the basis for the calculation of q(t) in Eq. 5. The analytic phase differences between channels yielded a symmetric, leptokurtic distribution (Fig. 5, dark curve), showing most of its values between Ϯ0.1 radians and long tails bounded by Ϯ radians. This differed from the nearly Gaussian distribution of phase differences derived from the Hilbert transform of band-pass filtered noise (dashed curve). The leptokurtic distribution reflected the low incidence of large changes in phase differences in the analytic signals from the EEGs. Figure 6 (top frame) shows histograms of the EEG values for the synchronization index, q(t) from Eq. 6, in the 4 solid curves, which show the means for the 4 cortices (visual, auditory, somatosensory, and entorhinal) pooled over all of 12 experiments on 2 cats. The percentages calculated in the following text refer to a mean of the four traces. The histograms reveal a high incidence of values above 0.5. The 99th percentile of the mean shuffled control data (dashed curves) is 0.59, and 28% of the area under the data histogram is above this value. The 99th percentile of the distribution for filtered noise (dotted curve) is just above 0.2.
The histograms of the spatial variation index, SD x (bottom frame), have a sharp peak at the low end of the distribution (solid curve), with 33% of the area below the 99th percentile (0.075) for the shuffled noise (dashed curve), and 70% of the area below 0.1. The cutoff for noise is 0.3. Comparing the two indexes of synchronization, the relatively sharp peak of the SD x histogram indicates that it may provide a better-defined index of intracortical synchronization. The bottom frame shows the distribution of the variation of the phase, SD t . This histogram is included to demonstrate that the variation of the individual channel phases is greater than that of the phase differences between pairs of channels. It shows that synchronization between channels is maintained while the individual channels' phases vary. It is not a trivial phase locking such as would be observed for two signals with almost the same constant frequency.
Another example of this result in 10-ms steps from a single trial for the visual cortex of a cat (Fig. 7, top frame, solid curve) shows that values of the normalized synchronization index, q(t), as a time series are consistently higher than the upper 99% confidence interval for the index computed from noise (bottom dashed line) and intermittently higher than the FIG. 6. Top frame: histograms of the distributions of q(t) from 4 data sets (solid curves) in comparison to values from random numbers (dotted curves) show the high values reflecting sustained synchrony. The shuffled data sets gave intermediate distributions (dashed curves). Middle frame: the spatial variation, SD x , gave a distribution that was clearly separated from the distributions of both random numbers and the shuffled data. Bottom frame: the distribution for the temporal variation, SD t , clearly differed from the distribution of random numbers but not from that of the shuffled data.  Fig. 3 (top frame). The curves show that the differences measured at 2 ms intervals do not change instantaneously to new levels but in the time interval of several msec. This type of display shows that the jumps cannot be attributed to problems with the method of unwrapping.
FIG. 5. A histogram is shown of the distribution of phase differences between channels, ⌬ x1,x2 (t), after normalization to remove the mean differences shown in Fig. 2. The sets for each 80-ms moving window were pooled to give the means (dark curve) Ϯ SDs (dotted curves) and to compare them with the distribution of normally distributed random numbers processed in the same way (dashed curves). The high central peak at zero in the range of Ϯ reflects the dominance of the plateaus in Fig. 3, top frame, giving the small values of ⌬ x1,x2 (t), and the broad distribution of the less frequent jumps predicted from Fig. 1. 99% confidence level for the index from the shuffled data (top dashed line). The index for these intracortical EEG data shows intermittent drops from a sustained high level, but not so the shuffled data. The decoherence index, SD x (t) (bottom frame, solid curve) shows peaks in the variance that are almost mirror q(t) in the top frame. The temporal SD, SD t (t), had peaks that often coincided with the peaks of SD x (t) but not always and then with lower extremes of variance (bottom frame, dashed curve). The two indexes, SD x of synchrony and q(t) of decoherence, give comparable although mirror-image results, in which q(t) emphasizes maintained levels of synchrony, while SD x (t) illuminates the episodic departures from synchrony.
The relations were investigated between the rates of recurrence of the peaks of SD x and the times of onset of phase cones in the same data sets. Figure 8 shows the interval histograms in which nine data sets from six rabbits were pooled after visual inspection of the individual distributions. The minimum duration of the interval for the phase cone data was set by the minimum length of the search window, which in turn was guided by optimization of classification of AM patterns based on use of phase cones for localization (Freeman and Barrie 2000). The distribution of the phase decoherence episodes was not bounded in this way. Figure 9 shows the power spectral densities of the autocovariances of the unfiltered EEGs, of the phase cone sequences represented by step functions to indicate the estimated times of starting and ending of each cone, and the time series of the synchronization and decoherence indexes from eight representative data sets. The autocovariances were from 3-s segments Ϫ1,000 to ϩ2,000 ms with respect to stimulus onset and averaged over 40 trials in one session. Peaks were found in and below the theta range for the four time series in some trials and subjects, not in others, and not at identical frequencies in the same sets. Figure 10 shows cospectra between the data sets, which indicate some degree of shared power across the theta range. It was apparent that the evanescence of EEG theta activity and the aperiodicity in the recurrence of the cones and the decoherence episodes precluded finding close statistical relationships among these variables with this linear technique.

D I S C U S S I O N
These results using the Hilbert method conform in several details to results predicted from Fourier-based methods on the same EEG data. Whereas Fourier decomposition with cosines (Freeman and Barrie 2000) or with frequency-modulated wavelets (Freeman and Baird 1987;Freeman and Viana Di Prisco 1986) starts with specification of a frequency at which to calculate the average phase over a time segment, Hilbert decomposition begins with each point in a time series to calculate an analytic amplitude, v, and a phase, , where the frequency is implicit in the rate of change in phase or is calculated as a mean value over a time segment. Each method has advantages and disadvantages. The Fourier methods are well-suited to detect spatial patterns of phase when the rate of change in phase is sufficiently small over a time segment that the frequency can be approximated by a single mean value, but they are ill-suited to resolve rapid changes in phase and/or frequency. The Hilbert method gives sharp temporal resolution, but it is sensitive to corruption by noise, and, on signals with broad spectra, it gives seemingly random variations. The two features of these EEG data that made it possible to apply the Hilbert transform were the prior findings that the spatial AM patterns relating to behavior were most clearly identified in the identified gamma ranges of rabbit and cat EEG data, which justified temporal band-pass filtering, and that stable spatial patterns of phase in the form of cones existed in the rabbit EEG data recorded simultaneously from 64 electrodes at high density on each cortex, which justified partitioning of the phase variance into temporal and spatial components, SD t and SD x . Top frame: the synchronization index, q(t), was calculated in 80-ms frames stepped at 10-ms intervals for all pairs of channels from the spatial phase differences, ⌬ x1,x2 (t), in each cortical data set as a normalized function ranging (Eq. 6) bounded between 0 (no synchrony) and 1 (complete synchrony). The values of q(t) derived from EEGs (solid curve) showed sustained levels of synchrony interrupted by dips toward the noise level (bottom dashed line) reflecting loss of synchrony in time periods surrounding the phase discontinuities (Fig. 3, bottom frame). The shuffled data (top dashed line) did not show these dips. Bottom frame: the time series are shown for the SDs of the spatial phase differences, SD x (t) (solid curve), and of the temporal phase differences, SD t (t) (dashed curve), for the same data as in the top frame. These indexes of variation gave curves that were close to mirror images of the synchronization index, q(t).
The analytic method gave increasing or decreasing phase differences with reference to the zero phase difference at any arbitrary starting time of the computation (Fig. 3, top frame). The sequential phase differences revealed a staircase pattern in time resembling a random walk process, with sudden changes demarcating plateaus that had decreased temporal phase variation and stabilized spatial phase differences. The distributions of phase differences between channels from whole trials were dominated by the values within each plateau (Fig. 5), giving the narrow peak between Ϯ0.1 radians (Ϯ18°). These low values were consistent with the gradients of phase cones measured within array windows. The derivative, d/dt, showed rapid changes that occurred on multiple channels, not simultaneously but comparable to the 4-ms segment across which the derivatives were computed numerically. This narrow distribution was compatible with the hypothesis of a brief time delay in the spread across each cortex of a state transition marking the onset of each new wave packet. The amplitudes of the step changes varied widely, as expected for numerical differencing close to the 2-ms digitizing step. The cumulative but not monotonically increasing differences in phase in the "staircase" pattern were consistent with the random variations of the signs and locations of the apices of the phase cones. The histograms of the intervals between maxima in the decoherence index were similar to the histograms of the intervals between onsets of phase cones. The power spectral densities and cospectra of the EEGs, phase cones, synchronization index, and decoherence index, intermittently showed peaks in and below the theta range but without clearly defined relationships in the frequency domain. This outcome was consistent with the aperiodic recurrence of the underlying state transitions, leaving unresolved the demonstration of the precise time relations between the onsets of phase cones and the spikes in the phase derivative. That demonstration is important to relate the jumps in analytic phase to behavioral measures, as a method for confirming that the staircase appearance should not be attributed to filtered noise and to a random walk.
Another concern was that sudden jumps might be attributed to problems with unwrapping. Display at high temporal resolution revealed that the jumps tended to occur over several FIG. 9. Power spectral densities (PSD; solid) and SDs (dashed) were averaged over a representative trial set from each of 8 rabbits for the spatial ensemble average EEG, the phase cones represented by square waves, the phase decoherence index and the synchronization index. The large SDs reflect the wide variation in height and frequency of the peaks in PSD, especially the EEG with high peaks in 3 subjects near 4 -6 Hz and little or none in the others.
FIG. 10. The coherences from the co-spectra indicated shared energy in and below the theta range but failed to confirm significant coupling among the 4 time series.
digitizing steps in the same direction and seldom as single large jumps. They did not occur simultaneously at multiple channels but temporally dispersed, which was consistent with phase cone gradients and with the volume conduction of the EEG field potentials around each recording site.
The interpretation is suggested that a spike in the derivative of the analytic phase, d/dt, when it occurs on multiple traces, demarcates a cortical state transition (Freeman 2000b) from a diastolic (receptive) state to a systolic (transmitting) state (Freeman 2000d), after which stable AM and PM patterns are maintained for 80 to 100 ms. The progress of the group wave across the cortex under the electrode array is manifested at each electrode by the potential differences created by neurons in its vicinity, for which the half-amplitude distance in accord with the spatial spectra of neocortical epipial EEGs is approximately that of the interelectrode distance, 0.6 mm. The rapidity of change despite the blurring by volume conduction suggests that the change at each site is by a first-order state transition, implying the occurrence of a temporal discontinuity in cortical dynamics. The reverse state transition to a diastolic state appears to be more gradual. The spatial SD of the analytic phase differences, SD x (t), is better suited for setting thresholds to detect episodes of decoherence and possible discontinuity, whereas the synchronization index, q(t), as well as linear methods such as the Fourier transform, are more useful for detecting plateaus following the state transitions demarcating phase cone onsets. The time course of the synchronization index in Fig. 7 for the data in comparison to that for the shuffled data suggests that neocortices tend to maintain high levels of internal synchrony that are relatively briefly interrupted at irregular intervals by decoherence.
The lack of strong correlation between the maxima in the SD x (t) and the occurrences of phase cone onsets highlights further limitations of both methods. In particular, although the jumps in phase values at the apparent discontinuities in temporal phase sequences are bounded by Ϯ as indicated in Fig.  5, the magnitudes of those jumps are predicted to be homogeneously distributed in that interval, so that many of them may have been below the criteria for detection and therefore been missed. So also some phase cones may escape detection, if the criteria for selection are too stringent, and even if they are detected, the precision of determining starting times is poor, owing to the coarse-graining of the temporal windows in which the search is made for phase cones. False-positive identifications may result in both methods from excessive noise and irregularities in the data of unspecified kinds that enhance phase slip. Further refinements in the analytic method applied to EEGs may yield better time markers for the self-organized systolic wave packets that may mediate perception, but at present their precise temporal localization remains elusive. As its usage becomes clarified, the Hilbert transform may prove superior to the Fourier and wavelet methods currently used for the detection of phase cones in EEG data, as well as for localization of AM patterns using the analytic amplitude.
The organization of neuron populations in local domains having a common wave form of phase modulation with minimal temporal variance, SD t (t), and a stable, delimiting spatial pattern of phase with fixed spatial variance, SD x (t), constitutes a novel solution to the "binding problem," which in its linear form (von der Malsburg 1983) holds that the pulse trains of selected feature detector neurons in sparse networks transiently synchronize their timing at zero lag (Singer and Gray 1995). In contrast, the synchronization manifested in the EEG is at zero lag commonly seen in time averages. Moreover, the neurons involved are not limited to those directly activated by predetermined stimuli, and the AM pattern that is expressed by the mesoscopic activity relates not to the immediate stimulus but to the context of goal-directed behavior in which that stimulus has been learned in repeated training trials (Freeman 1999). According to this nonlinear hypothesis, the cortex is destabilized by diastolic input that increases the interactions among cortical neurons by virtue of their nonlinear gains (Freeman 2000a), so that the form of its systolic output carried by cortical action potentials is dependent less on sensory input than on massive local feedback from neighboring neurons, and the content of that feedback is determined more by synaptic modification during past and present learning than by immediate input. Readout is by divergent-convergent pathways that implement a spatial integral transformation that extracts the cooperatively generated carrier wave with its AM and PM spatial patterns, while attenuating the spatially incoherent, sensory-driven residues as no longer of interest or value (Freeman 2000a). In both hypotheses the outputs of cortices are carried by pulse trains on axons of distributed length, size, and conduction velocities; the larger the phase variance imposed by variance in axonal conduction delays, the smaller must be the output. The necessity for axonal conduction delays between cortices and their targets of transmission is not easily made compatible with zero time lag synchronization (Freeman 2000c), so the study of cortical dynamics requires continuing development of effective techniques for measuring the phase distributions of cortical activity patterns, to provide the data needed to evaluate models that address the problem at either the microscopic cellular level (Schillen and König 1994;Traub et al. 1996;Usher et al. 1993) or at the mesoscopic population level (Freeman 2000d).