| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Innovative Methodology
1 Center for Neuroscience, University of California, Davis, California 95616; 2 Department of Neurology, University of California School of Medicine, Davis, California 95616
Submitted 2 September 2003; accepted in final form 27 November 2003
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
), the planning and execution of movement (Donoghue et al. 1998; Fetz 1997), spatial memory (O'Keefe and Recce 1993), and other largescale neuronal processes (see Engel et al. 2001; Varela et al. 2001 for reviews). However, this subject remains controversial (Gray 1999; Shadlen and Movshon 1999; Singer and Gray 1995).
Although many classes of neurons are able to fire rhythmically by virtue of their intrinsic biophysical properties (Hutcheon and Yarom 2000), neuronal oscillations in vivo are often large-population phenomena. Understanding the mechanisms and the functional role of oscillatory activity in the brain requires the analysis of multiple time series, often combining different types of recordings, such as spike trains from several microelectrodes, local field potentials (LFPs), scalp electroencephalograms (EEGs), magnetoencephalograms (MEGs), and electromyograms (EMGs). Insights into mechanisms and function can be gained by identifying the subsets of neurons involved in these oscillatory epochs as well as characterizing the interactions among these subsets or between subsets of neurons and muscles.
Statistically significant correlations are interpreted as implying the presence of interactions between the neural populations recorded. In some experimental situations, the statistical reliability of detection can be improved by collecting multitrial data that are time locked to external events. However, in many cases, the statistical benefits of multiple trials are not available because the interactions under scrutiny are not event related but are rather spontaneous "ongoing" events or occur on a single-trial basis. The situation is worsened when the interactions to be detected are of short duration. In such cases, the compromise between statistical power and temporal resolution imposes a severe limit to detectability. In the present study, we present a statistical method to reliably evaluate interactions in such "single trial" cases.
Two approaches have dominated the analysis of temporal correlations in neural data, focusing on either the time or frequency domains. When the data consist of spike trains obtained from a series of experimental trials, the most commonly used measure of timing correlation is the crosscorrelogram (Brody 1999b; Kirkwood 1979; Konig 1994). In a multitrial experimental design, statistical significance is usually evaluated using a shift predictor across trials. The crosscorrelogram can be particularly useful to evaluate the presence of fast synaptic interactions and to measure time delays between spiking cell pairs, although in some important cases, this method can lead to spurious detections. For example, artifactual sharp peaks might appear for spike trains if there are slow covariations in the cells' membrane potential (Brody 1998). A second approach, most commonly used with continuous data, is to analyze correlations in the frequency domain by using the coherence spectra. Frequency-domain methods have some advantages over time-domain methods, particularly when the data contain oscillatory components. In addition to providing a clearer characterization of the processes underlying oscillations, the statistical properties of several spectral measures are well known, at least for stationary processes. Moreover, statistical spectral methods can be extended to point processes and might be used to address the questions of synaptic connectivity and time-delay estimation (Jarvis and Mitra 2001; Mitra and Pesaran 1999; Rosenberg et al. 1989). In both approaches, the common idea is to obtain a correlation coefficient, either in the time domain (i.e., as a function of time lag between the processes) or in the frequency domain, and evaluate its statistical significance. However, standard methods of statistical evaluation in both cases require the data to be stationary during the analysis period, a requirement that is rarely satisfied by neural data. To get around this problem, it is usual to use sliding analysis windows of short duration so the correlation coefficient (or coherence) can be obtained as a function of time. Time-frequency methods have been developed where statistical confidence intervals can be imposed on the coherence estimates, and they have been applied to a variety of signals, including spike trains and LFPs (Mitra and Pesaran 1999; Pesaran et al. 2002).
The subject of this study is a set of analysis tools of more recent development that are aimed at detecting correlations in the phases of the oscillators while discarding the effect of amplitude fluctuations (Feige et al. 2000; Lachaux et al. 1999; Tass et al. 1998). This approach is especially useful when the signals contain a clear periodic or quasiperiodic component and has been, in part, inspired by recent developments in the field of nonlinear dynamics, where phase synchronization has been a subject of intensive research (see, for example, Schafer et al. 1999). A number of modeling studies in this field have shown that phase synchronization of coupled nonlinear oscillators might occur even in the absence of amplitude correlations, particularly in cases where interactions are weak (reviewed in Pikovsky et al. 2001). A sensible approach to synchrony detection in oscillatory neural signals should therefore take into consideration that timing correlations are better evaluated in the phase of the oscillators. In contrast, the analyses of synchronization based on the crosscorrelogram or the coherence spectrum detect both phase and amplitude correlations and can therefore miss the type of event where correlations exist in the timing alone, while the amplitudes vary independently. The opposite situation, where covariations in the amplitude envelope of two signals can lead to artifactual assessment of timing synchronization, has been described in cross-correlation analysis of spike trains (Brody 1998, 1999a).
In this paper, we address the problem of detecting phase synchronization in pairs of neural signals that are naturally restricted to a narrow frequency band. The methods presented rely on the definition of a phase variable from the time series. Correlation statistics are drawn from the phase evolutions instead of the original time series. Phase synchronization is evaluated over time by computing three different correlation indices using a moving analysis window. A statistical criterion based on surrogate data sets has been implemented to detect the presence of significant correlations. We present a novel approach for generating surrogate phase evolutions from the data and compare it to existing methods.
We applied this method to tremor time series from EMGs and recordings from the basal ganglia (the internal segment of the globus pallidus; GPi) of parkinsonian patients. Previous studies of tremor cell pairs and neural-muscle and muscle-muscle pairs using classical spectral methods, suggest that tremor is generated from oscillatory activity that propagates through parallel circuits spanning central motor networks (Hurtado et al. 1999, 2000). This conclusion was reached by showing that in either tremor cells or muscle pairs, activity is significantly coherent between some of the oscillatory cells or muscle pairs but not others. For example, a tremor-related cell in GPi could be coherent with a leg muscle but not arm muscle even though both limbs exhibited tremor. Two cells, >3 mm apart, could be undergoing simultaneous oscillation but not be coherent to each other. These observations are consistent with the existence of separate clusters of phase-locked cells associated with particular muscles, which can display independent oscillatory activity.
The methods presented here have allowed us to examine the dynamics of phase-locking process within the pallidum and between pallidal neurons and muscles. An outstanding question is how phase locking between the parallel circuits in the basal ganglia-thalamocortical networks occurs as a function of time as this can reveal the characteristics of functional connectivity and the interactions between segments at key locations in the network. A similar problem of detecting time-dependent phase locking can be applied to other areas of neuroscience research where rhythmic activity is prominent and where the relative phase between spiking cells and other ongoing neural signals (LFP, EMG) displays dynamic changes over time. A good example is the study of phase precession in the spiking of pyramidal cells relative to the hippocampal theta rhythm (Harris et al. 2002; O'Keefe and Recce 1993).
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
, 2001). For the phase construction to be meaningful, it is important that the time series is oscillatory. The presence of oscillatory activity can be evaluated prior to the phase reconstruction procedure by inspecting the power spectrum of the data for significant peaks. Because most neural signals are nonstationary, and might present several time-varying frequency modes, it is recommended that the power spectral calculations be performed over time, using a lag window and/or multitaper method (Jarvis and Mitra 2001; Mitra and Pesaran 1999). Statistically significant oscillations can then be obtained by applying a jackknife estimation method (Thomson and Chave 1991). Once the significant frequency modes of the signals and their time locations are known, the oscillatory components at the frequencies of interest can be extracted by band-pass filtering the signal at a narrow band centered at each frequency band. The result of this procedure is an oscillatory signal that can then be utilized for the phase extraction procedure. The narrow band filter can be applied to different types of signals. In the particular case of rhythmic spike trains, the band-passed oscillatory signal is equivalent to a smoothed version of the spike train restricted to the frequency band of interest.
From the oscillatory signal so obtained, an angle of rotation, or phase, can be conveniently defined in the complex unit circle. Here we have adopted the Gabor method for phase reconstruction based on the Hilbert transform of the data (Gabor 1946). Alternative phase reconstruction methods have been reviewed elsewhere (Le Van Quyen et al. 2001; Rosenblum et al. 2001). Phase reconstruction methods can be applied to a variety of signals; here we consider neuronal spike trains and EMG, but it can also be used to study phase correlations in MEG, LFP, and EEG and signals (Gross et al. 2000; Rodriguez et al. 1999).
The phase evolution variables are then used to compute correlation indices that assess the strength of phase locking between the two signals. We have tested three different types of indices: phase coherence (Feige et al. 2000; Lachaux et al. 1999), phase difference entropy (Tass et al. 1998), and mutual information (Palus and Hoyer 1998). The indices were normalized to return a value equal to 1 for fully phase-locked signals and zero for completely uncorrelated signals. Note that these indices are only meaningful if the phase reconstruction is valid, a condition that is satisfied when the signal contains a clear oscillatory component. We have therefore used a signal-to-noise ratio (SNR) parameter, computed over time, to establish the presence of oscillations and discarded from the analysis data segments where either of the signals in the pair falls below a criterion SNR.
To draw meaningful conclusions from the indices, we need to know the distribution of their values that would be expected if the two signals were independent. Only index values that are on the high end (e.g., 95th or 99th percentile) of the distribution can be considered to indicate the presence of phase locking. That distribution is unknown, and we rely on surrogate data to draw an approximation to it. Importantly, the shape of the distribution, the associated cutoff values, and, therefore, our conclusions will depend crucially on the assumptions that are made in generating the surrogates (Palus and Hoyer 1998). In the following text, we compare the usual approach of obtaining surrogate data from a Gaussian process with the same mean and SD as the data with a novel strategy for creating phase surrogates from the instantaneous frequency series. The outcomes from different surrogate data schemes are compared for a pair of independent oscillators and a pair of coherent ones, obtained from neurophysiological recordings.
A crucial parameter in the statistical analysis is the number of oscillation cycles that are considered in the calculation. Longer observation times (i.e., more oscillatory cycles) will return more reliable estimates provided that the time series remained stationary throughout the observation period. If the time of observation is too short, important interactions can be missed especially if they are weak or masked by noise. On the other hand, if the data are nonstationary during the recording period, the indices obtained from long observation times might fail to capture the presence of transient interactions. To more fully evaluate the data, we performed the analysis for several different observation lengths, i.e., multiple time windows of varying durations, ranging from 1.5 to 15 s, equivalent to 6–60 oscillation periods at the center frequency of interest (4 Hz).
The statistics are therefore calculated for each point tk in time, taking N prior points in the phase series (which constitutes a window for analysis) to obtain an estimate of the phase-locking statistic. The analysis is repeated for different values of N, to cover the range of analysis windows utilized. The resulting estimates provide more or less time accuracy and statistical reliability. The data-analysis methods were implemented in MATLAB. The finite impulse response (FIR) filter design and application routines were obtained from the MATLAB Signal Processing Library.
Data collection
We used data from four patients selected from our on-going studies of patients with idiopathic Parkinson's disease (PD) and prominent tremor in the extremities.
Electrophysiological recordings from the pallidum were obtained from patients A and B during pallidotomy for the treatment of the patient's parkinsonian motor symptoms. At the time of surgery, patients had been off antiparkinsonian medication for 
8 –10 h. Informed consent was obtained from all patients. The protocol for collection of intraoperative data were approved by the Institutional Review Boards of The University of California and Kaiser Permanente Research Foundation.
A Radionics stereotactic frame was used for target localization. The coordinates of the target site for the first electrode pass, the optic tract just ventral to the ventral border of GPi, were determined from magnetic resonance images of the brain in relation to fiducial markings on the stereotactic frame. High-resolution images were also used to determine the relative location of optic tract, pallidum, and internal capsule. A craniotomy (2 cm diam), 2 cm lateral to the midline and 2–3 cm anterior to the coronal suture, was performed under local anesthesia to provide an entry site for the recording electrodes. The microelectrode was passed through an outer, stainless steel guide tube. The electrode trajectory was approximately parallel to the sagital plane at an angle of 30° from vertical. The initial recording track was set as 0,0 and the target site initially set at a depth of 75 mm from the surface of the cortex. Recording began at 2 cm below the cortical surface, and the electrode was moved slowly along the track using a micromanipulator. The purpose of the microelectrode recordings was to define the location for lesioning (the ventroposterior part of internal segment of the GPi) and to identify the safe borders for lesioning (avoiding the internal capsule and optic tract). According to the Schaltenbrand and Wahren (1977) atlas, at a laterality of 21 mm, GPi is at 64–65 mm and optic tract at 75 mm.
Microelectrode recordings were obtained with 50-µm, beveled, stainless steel electrodes (0.5–1.0 M
impedance at 1 kHz). Signals were amplified (10,000x) and band-pass filtered between 500 Hz and 10 kHz. EMGs were recorded from two arm muscles [abductor pollicis brevis (APB) and wrist extensor] and two leg muscles (quadriceps or tibialis anterior and gastrocnemius) contralateral to the side of surgery with Grass disk electrodes, amplified (10,000x) band-pass filtered at 30 Hz to 1 kHz. Electrical signals from the recording electrode and EMGs were stored on analog tape using a Vetter recording adaptor and VCR. Neuronal and EMG signals were digitized using an RC Electronics (Santa Barbara, CA) A/D board on an IBM-compatible, Pentium-based, personal computer running RC Electronics software. EMG and brain signals were digitized off-line at 10 kHz. Spikes were threshold extracted (>2 SDs above baseline). After spike extraction, spike trains were down sampled at 1 kHz.
Data are also presented from two patients (C and D) with prominent resting tremor who were among a group of PD patients studied using EMG recordings only. Informed consent was obtained from all patients. The protocol for collection of EMG data from patients C and D was approved by the Institutional Review Board at The Veterans Administration's Northern California System of Clinics.
In this group, we recorded from up to eight muscles. EMG signals were obtained with disposable silver/silver chloride electrodes (Multi Bio Sensors). Drug therapy was not withdrawn prior to the EMG recordings. Electrodes were placed in bipolar arrangement, each pair positioned along the axis of the recorded muscle, with 4 cm separation between leads. EMG signals were differentially amplified 1,000x with a Myosystem 2000 (Noraxon, Scottsdale, AZ) system, filtered (16–500 Hz), and digitized at 500 Hz with a 16-bit data acquisition software (Datapac, Run Technologies, Laguna Hills, CA) running in a PC.
Instantaneous phase reconstruction from time series
To obtain a phase evolution variable from a time series, we first filter the signal in the band of interest (Fig. 1A). This reduces the effect of noise in the phase estimate. In the case of parkinsonian tremor, where the power spectrum has a typical peak at 3–5 Hz, we used a 512-pole band-pass FIR filter centered at 4 ± 2 Hz. Because we are interested in extracting phase correlations, phase distortions due to filtering must be corrected. To eliminate phase distortions, the data were filtered recursively in two steps; the second step uses the same filter parameters but is run backward in time. The filter design and implementation steps were performed with routines from the MATLAB Signal Processing Library.
|
After filtering, we utilize the complex analytic representation method introduced by Gabor (Boashash 1992; Gabor 1946; Schafer et al. 1999) to extract the phase evolution (Fig. 1B). The complex analytic extension of the filtered signal x(t) is given by
 | (1) |
(t) is given by the Hilbert transform of the signal
 | (2) |
(t) with the same power content as the original x(t) but phase-delayed at each frequency by –
/2 [e.g., the Hilbert transform of cos(
t) is sin(t), and its Gabor representation is the unit circle trajectory cos(t) + i sin(t)]. Instead of evaluating the integral in Eq. 2, it is more practical to obtain the Gabor representation in the frequency domain by taking advantage of the property that the spectrum of an analytic signal is identically zero for all negative frequencies. The analytic signal in Eq. 1 can then be easily computed by taking the discrete Fourier transform (DFT) of the data, setting negative frequency values to zero, and performing the inverse Fourier transform.
To obtain a phase evolution of the oscillation, the analytic signal is projected into the unit circle (Fig. 1B)
 | (3) |
(t)|| is the modulus of (t) and the phase 
obtained as the argument of the resulting complex series. For discrete data, the procedure yields a series of complex values of unit length zk = z(tk), at time points tk. The phase series is then obtained by taking the complex argument (angle) of the zk
 | (4) |
It is convenient to define phase as the "unwrapped" angle: instead of being periodic, the phase advance accumulates over time so that a cycle completion is 2, two cycles 4 and so on (Fig. 1C, top). From the phase so obtained an angular frequency can be defined as
 | (5) |
t = tk – tk–1 is the sampling period. We exploit this relationship later to generate surrogate phase series by integrating surrogate instantaneous frequency series
 | (6) |
A better estimate of the instantaneous frequency can be obtained by applying a higher order differential filter or, equivalently, by smoothing the preceding frequency estimate (Boashash 1992; Schafer et al. 1999). The meaning of an instantaneous frequency and its relation to the spectrum of a signal as well as numerical methods for its computation is treated extensively in (Boashash 1992). Although the expression in Eq. 5 is a first-order approximation that can give rise to nonsmooth instantaneous frequency traces, we are not interested in estimating instantaneous frequency per se but in analyzing the coupling of paired recordings. We utilize the instantaneous frequency as an aid to extract phase slip events from the data and to construct surrogate phase series that are then used to evaluate the statistical significance of phase-locking indices, procedures that are described in a later section. This does not require higher-order estimates. Indeed, as will be apparent later, the construction of surrogates relies on time integrating the instantaneous frequency.
Phase singularities: aperiodic data, phase slips
PHASE SINGULARITIES. When applied to quasisinusoidal signals, the Gabor method yields a continuous angular rotation variable. Its first derivative, the instantaneous angular frequency, is also continuous and should be approximately bounded within the values of the frequency band of the input signal, which is determined by the band-pass filter parameters (Boashash 1992; Rosenblum et al. 2001; Schafer et al. 1999). However, in experimental data, periodic activity is often interrupted by episodes of aperiodic activity (Fig. 2) or by phase slips (Fig. 3), resulting in discontinuities in the phase evolution variable. Interruptions can be brief, as is the case with phase slips (Fig. 3), a situation that we treat in detail in the next subsection, or of longer duration, when oscillatory activity ceases altogether (Fig. 2). In both cases, phase discontinuities appear as "spikes" in the instantaneous angular frequency plot. The particular shape and amplitude of the frequency spikes does not reflect any intrinsic property of the oscillatory series because they depend, as well, on the particular filtering parameters and method utilized to compute phase.
|
|
PHASE SINGULARITIES ARISING FROM APERIODIC DATA. Angular frequency spikes occur when a signal ceases to be oscillatory and becomes noisy or aperiodic. The power spectrum of such signals is characterized by a broad spectral spread. In these cases, the definition of phase loses its meaning. Because detection of phase locking requires a reliable phase construction, we must determine whether the signal is oscillatory before proceeding with the analysis of phase correlations. Intervals where one or both signals are not oscillatory are not considered in the analysis of phase locking described here.
As a criterion to discriminate between oscillatory and aperiodic segments we define the signal to noise ratio (SNR) as the ratio of power in the band of interest over the total signal power
 | (7) |
, t) is the "instantaneous" power at frequency and a, b are the limits of the tremor band (2–6 Hz). P(, t) is obtained by computing the DFT of the signal using a sliding window of 1-s duration. To improve reliability, P(, t) is taken as the average of 10 successive windows with 90% overlap. By setting a threshold value for the SNR, we can pull out the segments of well-behaved oscillations and perform the phase-locking analysis on them. On the basis of inspection of many tremor time series, we found that a threshold a value of SNR >3.7 is a good indicator of oscillatory behavior for our data. We therefore used in the analysis only data segments satisfying this criterion. Segments with lower SNR were discarded as their phase representation is not meaningful. A comparison between the phase representations of an oscillatory signal and an aperiodic one is shown in Fig. 2. Figure 2A shows a time series of a tremor signal recorded from an arm muscle, which is interrupted by a short episode of "noisy" or aperiodic myoelectric activity. During the oscillatory episodes, the Gabor representation of the signal is a smooth rotation, of varying amplitude, around the origin. When these trajectories are projected onto the unit circle, their angular frequency is a continuous function of time. In contrast, during the episode of aperiodic activity, the trajectory is no longer rotational but crosses the origin, resulting in "jumps" in the unit circle projection (Fig. 2B). The derivative of the phase gives rise to singularities or spikes in the instantaneous frequency plot. Although the filtered signals of both the noisy (gray line in Fig. 2) and oscillatory segments (black line) are restricted to a narrow frequency band (4 ± 2 Hz), only the oscillatory component has a well-defined phase.
PHASE SLIPS. Phase slips are singularities or jumps in the phase evolution of an oscillator. They may arise, for example, in forced oscillators when the forcing function frequency differs from natural frequency of oscillators or in networks of coupled oscillators with mismatched frequencies. Phase slips have been extensively described in a number of modeling studies and in experimental studies, e.g., in phase resetting phenomena (for examples, see Glass and Mackey 1988). Here we present a simple method to detect phase slips in experimental time series and to quantify the amount of phase advance associated with particular events.
Phase singularities are not artifacts of the particular method used to compute the phase advance nor are they a result of the particular band-pass filter used. Although all these factors do affect the particular shape of a frequency spike, in studies of phase locking, we are only concerned with the time location and the amount of phase advance associated with the occurrence of the phase slip, which should be independent of the particular phase reconstruction method utilized. Note that at a singularity, the instantaneous angular frequency reaches values that are well above or below the limits imposed by the band-pass filter. This is the consequence of phases "jumping" as the complex trajectory crosses the origin. In contrast, phase is well defined when its derivative, the instantaneous frequency, is a smooth function of time.
To find the time location and amount of phase advance associated with a phase slip, we use the instantaneous angular frequency series computed as in Eq. 5, where phase slips appear as spikes. By setting as thresholds the band-pass limits of the filter, the times of slip occurrence can be found at the points of threshold crossing. The procedure is exemplified in Fig. 3.
The phase advance caused by a slip can then be quantified by comparing the phase immediately after the singularity to the forecasted phase if the singularity had not taken place. Taking advantage of the fact that the phase is the time integral of the instantaneous frequency, a forecasted phase can be defined by time integration of a version of the instantaneous angular frequency where the singularity has been erased. We assemble a continuous instantaneous angular frequency 
from the original series by erasing the singularity and filling the vacant points by cubic interpolation (Fig. 3A). is then integrated to construct the forecasted phase series 
(Fig. 3B)
 | (8) |
The phase advance due to the singularity is computed as
 | (9) |
Phase slips occur frequently in tremor-related activity in the basal ganglia and, therefore it is important to include them in the construction of surrogate data to estimate the distribution of phase-locking statistics.
Phase-locking indices
In this section, we consider the phase synchrony between a pair of oscillatory signals. The general approach to the detection of phase synchrony relies on computing a correlation statistic from the data. The goal is to test whether the relative phase, or phase difference, between the signals is bounded
 | (10) |
. This is a particular case of the more general n:m phase-locking case, where
 | (11) |
; Tass et al. 1998). Because we are studying signal pairs that have prominent spectral peaks at a similar frequency, we focus here on the 1:1 locking condition, but the methods can be modified to study other locking modes. All three phase-locking statistics considered in the following text have been normalized between 0 and 1. Zero corresponds to a pair of independent signals and 1 to perfect locking. Even though all statistics are normalized to 1, comparison of values between 0 and 1 should be made with caution because these three measures are not necessarily equivalent to each other.
PHASE COHERENCE. The phase coherence index is based on the [magnitude squared] coherence spectral estimator for bivariate time series (Rosenberg et al. 1989). We define the relative phase series as 
j 
1(tj) – 2(tj), j = 1... N, where the tk are the sampling points (see Fig. 4A). The time-dependent phase coherence is then defined as (Fig. 4B)
 | (12) |
, where fs is the sampling rate of the signals. The phase coherence is always less than or equal to 1, taking a value of 1 only when the relative phase remains constant throughout the observation period T. A similar measure has been used by others (Lachaux et al. 1999; Rosenblum et al. 2001).
|
; Tass et al. 1998) is also computed from the series j for each time point tk. A histogram of j for j = k–N,..., k the observation time window is first built. The entropy of the series is then defined as
 | (13) |
This value can be normalized to the maximum entropy, which is achieved for the uniform distribution: pj = 1/L for all j, and hN reaches its maximum value 
. The normalized entropy is then
 | (14) |
MUTUAL INFORMATION. A mutual information index is defined as
 | (15) |
, 
for l = k–N,..., k and pi,j is obtained from the joint histogram of the pairs , (see Palus and Hoyer 1998 for overview). This statistic is computed over the histogram of N accumulated phase values at each point in time. The procedure is illustrated in Fig. 4D where the joint histogram is shown in the grayscale density map.
The mutual information can also be normalized to its maximum value, IN = log L, achieved when the series , are identical,
 | (16) |
The mutual information has the additional advantage that it achieves high values for any n:m phase-locking mode, in addition to the 1:1 case. However, because we are prefiltering the data in the same range (tremor range), only 1:1 interactions will be expressed in the information index. In addition, mutual information picks up so-called nonlinear correlations between phases.
Note that the pairs , considered in the histogram are paired in time so that mutual information as computed here does not consider time delays between the time series. This is in contrast to other mutual information measures where time delay is set as a parameter (Hoyer et al. 2002). However, for the cases under study, where the data have periodicities, it is not crucial to consider the time delays because any delayed correlation will be expressed as "stripes" in the joint probability distribution.
Surrogate data schemes
To use the phase-locking indices in a meaningful way, we need to know their distribution under the null hypothesis of independent pairs of oscillatory activity. Only values that depart significantly from what would be expected for independent oscillators can be considered as revealing the presence of phase locking. In this paper, we are dealing with single-trial data and cannot shuffle trials to determine the null distribution. And because we have no prior knowledge of the series statistics, it is not possible to infer the distributions analytically. A way around these problems is to generate, for each of the data channels, a large ensemble of surrogate time series that share some statistical features with the original data (Palus and Hoyer 1998; Schreiber and Schmitz 2000). This approach has been applied successfully in several studies where a statistical assessment of a property is necessary but the number of observations are limited. The distribution of the index, computed for pairs drawn randomly from the surrogate ensembles, can be considered as an approximation of the null distribution. This approach has become more feasible in recent years with the increase in low-cost computing power.
The null distribution so obtained will depend crucially on the particular method used to generate the surrogates. Methods differ in the specific features of the original data that are kept fixed across the surrogate ensemble, say, mean and SD, or power spectrum. As more features of the data are preserved, the resemblance between the surrogate and the original increases, and the test is a stricter one.
Because the phase-locking indices are computed using sliding analysis windows of variable length, it was necessary to obtain a separate cutoff value for each one of the window lengths utilized. This is because the distribution of phase-locking indices for a pair of independent processes depends on the length of the data stretches entered into the calculation. If two independent signals oscillate at similar frequencies, their phase lag will remain nearly constant during one oscillation period, thus giving rise to a high phase-locking index value. As the observation time increases, and more oscillation cycles elapse, the phase differences will diffuse away from the initial angle, eventually covering all possible angular values. This would result in low phase-locking index values. In the ideal situation of infinite observation time, the phase-locking index for independent processes should be equal to zero.
We generated surrogate data at the length of the longest analysis window. The index values for different window lengths were calculated by appropriately increasing the span of the surrogate segment included in the calculation. This avoided the need for generating a new set of surrogate data for each window length utilized.
We tested four different surrogate data methods. In the first method (S1), a surrogate time series is assembled by drawing samples from a Gaussian distribution that preserves the mean and SD of the original series. From the surrogate time series, the phase evolution and the coupling indices are computed as explained in the preceding text.
In the remaining three methods (S2–S4), we generate surrogate phase series directly instead of obtaining time series first. Constructing surrogate phases directly bypasses the steps of filtering and phase reconstruction, thus reducing the amount of computation. Because the phase-locking indices that we are testing are computed from the phase series, for our purposes we only need the phase representation of the surrogate data to resemble the phase representation of the original data. The methods we present here consist of generating an instantaneous frequency surrogate and obtaining from it the phase series by integration. We compare the performance of the different schemes by applying them to pairs of time series which we know to be independent. We created these by pairing the recordings of tremor activity from two different patients. For this type of data, any phase locking detected is a false positive. In addition we compared the methods by applying them to pairs of signals with high phase coherence averaged over a long time window.
SURROGATE METHOD S1: LINEAR GAUSSIAN AND POISSON PROCESS. For spike trains, random spikes are generated from a Poisson process with the same mean firing rate as the data. For other data types (EMG, LFP), surrogate time series are generated from a Gaussian distribution with the same mean and SD as the original data (Fig. 5A).
|
is given by 
, where t is the time resolution of the phase series and 
a step size that can be considered as a surrogate instantaneous frequency variable. The phase evolution is then constructed as
 | (17) |
of surrogate instantaneous frequencies. The approach has been to obtain the from the instantaneous angular series computed from the test data (Fig. 5, C and D). Here we present three alternative methods for selecting the .
SURROGATE METHOD S2: SURROGATES FROM RANDOM PERMUTATION OF INSTANTANEOUS FREQUENCIES. In this method, the surrogate instantaneous frequency, , is taken by random permutation of the original instantaneous frequency values (Fig. 5D, top), k. The rationale is that the surrogate phase evolution will phase-advance, on average, as much as the original data, but any phase correlation in the two series will be destroyed in the randomization process. Prior to the randomization process we removed the spikes (i.e., phase slips) in the instantaneous frequency series. As we stated previously, the shape of the spike that presents a phase slip in the instantaneous frequency series depends on the numerical method used to compute phase. On average resampling the spikes artificially increases the speed of phase advance.
SURROGATE METHOD S3: SURROGATES THAT PRESERVE THE POWER SPECTRUM OF THE INSTANTANEOUS FREQUENCY. This method attempts to preserve some of the temporal structure in the instantaneous frequency series. To do this, we create surrogates having the same power spectrum as the k series (Fig. 5D, middle). The method consists of computing the DFT of the k series, randomizing the phase spectrum while preserving the amplitude spectrum (with the restriction that negative and positive phase values remain symmetrical), and then computing the inverse DFT to obtain the . As in the previous method, we eliminate the spikes caused by phase slips prior to the computation of surrogates. The reason to do this is that the spikes have power at all frequencies. A better way to take into consideration the phase slips in the construction of surrogates is presented in the following scheme.
SURROGATE METHOD S4: AS IN METHOD S3 BUT INSERTING PHASE SLIPS. A problem with methods S2 and S3 is that phase slips have been omitted from the surrogate phase walk. Phase slips are important because they can promote phase locking, by advancing the phase of a slower oscillation, or can have the opposite effect, preventing two oscillators of similar frequencies from becoming phase locked. In this method, we inserted the phase slips at random locations into the surrogate time series (Fig. 5D, bottom) by using the following procedure: 1) extract the phase slips from the instantaneous frequency series as described in Phase singularities: aperiodic data, phase slips (Eqs. 7 and 8). Obtain the phase advance associated with each slip and create the set 
p, p = 1,..., M where M is the total number of phase slips. 2) Generate surrogate instantaneous frequencies from the "spike-free" k, preserving its power spectrum as described in the previous section. 3) Select surrogate phase slip times tp from a Poisson process with rate 
= M/T where M is the number of phase slips and T is the time duration of the series. 4) Randomly assign a value p from the set of phase slips to each slip time, tp. And 5) obtain the phase values from the using 
for k corresponding to tp and 
elsewhere.
Comparison of surrogate data methods
In this section, we compare the surrogate data schemes presented in the preceding text. To be accurate, the statistical method must minimize detection errors: false positives (rejecting the null hypothesis when it is true) and false negatives (accepting the null hypothesis when it is false).
To test how likely the different methods are to produce false positives, we computed the phase-locking indices on a pair of known independent processes, generated by pairing tremor series from two different patients. If a particular method yields a positive answer (i.e., it "detects" episodes of phase locking), we know by construction that it was a chance occurrence.
A comparison of the three surrogate data schemes, as applied to the independent tremor series, are shown in Fig. 6. The left column shows the 99% cutoff value of the three phase-locking indices (phase coherence, entropy, and mutual information) obtained from the surrogate data ensembles as a function of the integration time window. Note that the 99% cutoff value for each of the indices is high (i.e., near 1) for short integration times and approaches zero as the duration of the analysis window gets longer. This is because even for independent processes, if oscillations are of similar frequencies, the phase difference cannot change much over times that are short relative to the duration of a cycle. As more oscillation cycles are included in the analysis window, the phase difference will tend to cover a wider range of angles, resulting in lower phase-locking index values.
|
1% false positives. Yet the Gaussian surrogate scheme (S1) yields a much higher rate of significant values. This is also the case for method S2, which is similar to S1 in that it also destroys the temporal structure of the data, although it does so in the phase domain. In comparison, methods that preserve the power spectrum of the k (methods S3 and S4) did not return false positives for any of the analysis window lengths.
We are also interested in evaluating how likely the different methods are to produce false negatives, whereby interactions that are present are not detected because the cutoff criterion is too strict. A rigorous test for this would require data where the phase correlations and their time locations are known (i.e., by construction). However, we do not have a priori information on whether recorded pairs are or are not correlated. Therefore we selected a paired recording that has high coherence that is statistically significant according to a commonly used method (Hurtado et al. 1999; Miller and Sigvardt 1998; Rosenberg et al. 1989). The data from this example consist of a simultaneous microelectrode recording from the globus pallidus and an EMG from the medial gastrocnemius during a tremor episode. Figure 7A shows that both signals have sharp spectral peaks and are significantly coherent at the frequency of tremor, indicating a strong level of phase locking during the recording period. This phase locking can be examined over time using the phase-locking indices that we have presented above. In Fig. 7B, we show the time evolution of the three phase-locking indices, using two different sliding window lengths. A long integration window of 12 s (gray traces) yields index values above the 95% cutoff limit (gray dashed line) for most of the recording time. Note that each point in the curve shows the index value computed from the previous 12 s of data. The time locations where the index values are significantly high are marked by a gray rectangle. Note (Fig. 7B) that the rectangle precedes the points where the index trace crosses the 95% cutoff line. This is because the rectangle marks all the data points that contributed to the index values that exceeded the cutoff level.
|
More precise information about the time when phase locking occurs within the recording epoch can be gained by using sliding windows of shorter duration. The black traces show the time evolution of the same phase-locking indices but using a shorter sliding window duration (1.5 s), and the black dashed line marks the 95% cutoff value for that integration window length. The time segments where the signals are significantly phase locked are marked by the black rectangles. Note that by using an analysis window of shorter duration we have gained information as to which time locations show stronger phase locking (compare black and gray traces and rectangles). A more complete description of the characteristics of phase locking can be obtained by performing the analysis for a number of sliding window lengths. The periods of significant phase locking for each analysis window length are illustrated in Fig. 7C.
The three indices give similar answers for the S3 method and nearly identical answers when phase slips are included in the surrogate data (S4 method). Note that the time evolutions of the three indices display peaks at the same time locations, so they capture essentially the same features of the phase correlations. The differences in the time locations of significant correlations (Fig. 7C) occur because the peaks are close to the cutoff levels, and therefore significant points are sensitive to small changes in the particular selection of a percentile cutoff.
We have shown in the examples in Fig. 7 that the surrogate data method provides an effective way of detecting phase locking while preventing false positives (Fig. 6). In situations where the data series are largely stationary, the method we have shown is consistent with the standard coherence method, yielding significant phase-locking episodes that are sustained, or evenly distributed, over time if the series are coherent.
The advantages of this method become apparent when the data are nonstationary. An example of this is shown in Fig. 8. The recording consists of 30 s of simultaneous activity from a GPi spike train and a forearm flexor EMG during an episode of tremor. Both signals have sharp spectral peaks at 4 Hz (Fig. 8A), yet their coherence spectrum shows no sign of correlation at that frequency (the small peak occurs at 2 not 4 Hz). This result leads one to conclude that the signals are independent oscillators. However, the sliding window analysis (Fig. 8B) indicates the presence of phase-locking episodes. Shorter analysis windows (darker rectangles) pinpoint the time locations of highest phase locking (t1 and t2). Inspection of the distribution of phase difference angles during those two periods (Fig. 8C) shows that there was a shift (90°) in the phase difference between the two episodes. An analysis window long enough to include both episodes would yield low phase-locking values, giving the impression of no phase locking throughout that period. Indeed that is the case for the coherence estimate in Fig. 8A, where the full 30 s of data has been entered into the calculation. As the analysis window length is reduced, information is gained about phase-locking episodes spanning shorter time lapses. Using multiple time windows helps us to get around the problem of not knowing at which length scales the data remains stationary.
|
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
Here we have introduced a statistical method to detect phase correlations in single trials of paired neural data. The method applies to signals that are oscillatory in a narrow frequency band. To obtain a phase representation, we used a Gabor method, and performed the statistical analysis on the phase variables. We have restricted our analysis to segments that are oscillatory, like the ones encountered in tremor episodes, and have discarded episodes of irregular activity. The SNR measure we applied here serves the purpose of discarding activity that is broadband, where no oscillatory behavior is visible and where the phase representation is marked by discontinuities.
Phase reconstruction for spike trains
Because the phase evolution is defined as a continuous rotation in the complex plane, the question may arise as to whether or not it is legitimate to use the Gabor, or any other method, to construct a phase evolution for a point process such as a spike train. In this respect, it is important to note that the phase-reconstruction method described here is applied not to the point process itself but to a continuous signal obtained by running the point process through a band-pass filter. The filtered signal so obtained is akin to a spike density function (or instantaneous firing rate) that it is band-limited around the modal frequency of the oscillatory process. In fact, in most neurophysiological studies, the instantaneous firing rate is typically obtained by convolving the spike train with a windowing function (square, Gaussian, or other). Any such window will have a frequency domain representation with band-pass properties. Computing a spike-density function is therefore roughly equivalent to band-passing the original spike train. The FIR filter design used in this study is a particular choice of windowing that captures the signal component around the frequency of tremor activity, which is the dominant mode in the spectrum.
Such filtering does entail a loss of information on the timing of the individual spikes. What remains is a continuous signal whose peaks times can be interpreted as the locations of maximum spike probability within any given cycle. In the particular case where there is one spike per cycle, the peaks of the filtered spike trains should coincide with the time location of the spikes. Indeed, in some physiological point processes where there is one event per cycle, such as the electrocardiogram (EKG), the phase reconstruction can therefore be simplified without need of filtering, by assigning a 2 rotation between event times (see Schafer et al. 1999). In the case of trains of action potentials from neurons, there is no marking event per cycle and filtering is necessary. One may ask whether the loss of timing information of the individual spikes caused by filtering is acceptable. Indeed it is because the phase extraction is done at a particular frequency under the assumption that that particular frequency band carries the relevant information. In the case of a tremor time series, where signals are localized in the 2- to 6-Hz band, the assumption is a reasonable one. If there is more than one frequency mode, the analysis steps can be repeated for each of them by restricting the band-pass filter to the different frequency bands. In addition, phase coupling can be evaluated across different frequencies, but in that case, the phase difference has to be redefined as m, n = n2(t) – m1(t) where the center frequencies under consideration satisfy the relationshipnf2 = mf1 (Rosenblum et al. 1996; Schafer et al. 1999; Tass et al. 1998); the remaining procedure is the same as described in the preceding text. In summary, the method can be applied in succession to each the frequency bands of interest to extract all the relevant information in the spike train.
When applied to spike trains, the phase-reconstruction method can also be utilized to estimate the relative phase between a spiking neuron and a continuous signal such as the EMG, local field potential, or other. The instantaneous phase difference (t) = EMG(t) – S(t) provides an estimate of the (mean) phase lag of the spikes relative to the EMG. In Fig. 8C, we show histograms of (t accumulated over short time periods that reveal that the phase relationship can vary significantly across locking episodes. It might be questioned whether it is legitimate to obtain that phase relationship by using a filtered version of the spike train instead of the intact spike times. We have compared the mean relative spike phase obtained as in the preceding text with that obtained by Csicsvari et al. 2003. In that study, the phase of spikes relative to the ongoing theta rhythm is computed by assigning to each spike its phase value within the theta cycle and building an angular histogram from the cumulative spike phases. We used that method to compute the phase relationship between spikes and (filtered) EMG, and compared the mean of the histograms so obtained to the mean phase difference value obtained from our method. The values were essentially the same (data not shown).
PHASE SYNCHRONIZATION. Methodologically, the analysis of phase locking is a sensible way of detecting synchronization because it captures the covariations in the relative timing of the oscillations while neglecting amplitude fluctuations (Schafer et al. 1999). The importance of separating the analysis of phase and amplitude covariations in oscillatory networks can be better appreciated if one considers the results of recent modeling studies in the field of nonlinear dynamics. It has been shown in a variety of nonlinear models that as the strength of coupling between oscillators is increased, phase entrainment precedes the appearance of amplitude covariations (reviewed in Boccaletti et al. 2002; Pikovsky et al. 2001). For weakly coupled neural oscillators, phase entrainment might be strong even when the amplitudes are uncorrelated. A consequence of this is that a correlation index sensitive to both phase and amplitude covariations can be less effective in detecting phase locking than purely phase indices.
The method we implemented relies on phase-locking indices similar to ones that have been used previously but with two important innovations. First, we developed and tested novel surrogate data schemes that minimize the rate of false negatives, while retaining the ability to detect subtle correlations in the data. Second, we perform the analysis over a number of sliding window lengths, thus addressing both the need for large stretches of data to detect weak correlations and the need for short stretches to mitigate the adverse effect of nonstationarities. The cost of computing phase-locking indices for each of the surrogates is paid only once because one ensemble of surrogate data are used to find cutoff curves as a function of window length (Fig. 6B).
Surrogate data schemes
Surrogate data schemes were used to test the null hypothesis of independent activity. We introduced a novel approach whereby surrogates were assembled from a sequence of elementary steps around the unit circle, obtained from the instantaneous frequency of the original data. The advantage of this approach over methods that generate surrogates of the original time series is that it bypasses the filtering and phase extraction analysis steps, thereby simplifying the computations. In the simplest method (S2), the elementary steps were obtained from a shuffled version of the instantaneous frequency. This destroys the phase correlations between the series, but it also destroys the temporal structure of the phase evolution. We found that this approach yields a high level of false positives. A similar problem was encountered with the Gaussian (S1) method. In comparison, the methods that preserve some of the temporal correlations in the instantaneous frequency series (methods S3 and S4) provided a more realistic approximation of the dynamics and yielded a stricter cutoff level, eliminating false positives.
How strict a surrogate method is will depend on the expected time for the relative (surrogate) phase = 1 – 2 to cover the circle uniformly. For oscillatory signals with similar frequencies, the method that preserves the temporal structure of each signal yields a drift speed for that approaches more closely what would be the case for two real signals. In contrast, method S2 accelerates the drift of and gives a lower cutoff value for the phase-locking indices. It is important to take into account that the drift speed of depends largely on the particular data set under consideration. The most important factor is the difference between the center frequencies of the two oscillators: the larger that difference is, the faster the relative phase drift will be (
). Therefore shorter observation times are required to assess statistical independence of signals with disparate frequencies than for signals with similar frequencies. A second important factor is the spectral spread of the signals. For two signals with a common frequency mode, the diffusion of will be faster if the variance of the instantaneous frequency (i.e., spectral spread) is larger. In contrast, if the variance is small, the time needed to determine whether signals are phase locked or not will be large. A third factor is the temporal structure of the phase evolution. The phase evolution of the signals we have worked with has a smooth derivative, punctuated by occasional singularities (phase slips). Usually we find that the instantaneous frequency contains different trends over time. A realistic surrogate data construction should attempt to replicate these characteristics. We found that S3 and S4 surrogate methods, which preserve temporal structure, provide a satisfactory way of representing these trends.
In contrast, we found it convenient to extract fast events (i.e., phase slips) and treat them separately, either neglecting them (method S3) or adding them to the surrogate series as a poisson point process (method S4). If the phase slips are not extracted, the S3 method will tend to spread over time the power that is present in the singularities, resulting in a noisier series without slips. Including phase slips in the surrogate scheme has the effect of accelerating the diffusion of , thereby reducing the cutoff levels. Of course the effect will depend on the rate of phase slips in the data and in cases where slips are rare the differences between methods S3 and S4 are negligible.
Multiple analysis windows
The classic method to compute coherence subdivides the data into a number of nonoverlapping windows of fixed duration and utilizes the DFT of each window as a data point (reviewed in Carter 1987). That method has the advantage that statistical cutoff limits can be easily computed, using a formula that does not depend of the specific properties of the individual time series (Miller and Sigvardt 1998; Rosenberg et al. 1989). However, it assumes that the data are stationary during the observation epoch. Statistical significance can be gained by increasing the number of independent windows entered into the coherence estimate, but this requires either increasing the total time length of the data to be analyzed or reducing the duration of each window. Although statistical reliability can be increased by using overlapping windows (Carter 1987), stationarity is still required and the tradeoff between time resolution and statistical reliability persists.
These methods are problematic for most real neural data, which is nonstationary by nature. To detect transient synchronizations in nonstationary data requires short analysis windows, yet the statistics will only be reliable if the correlations are strong enough; weak correlations require longer observation times to be detected reliably. The approach presented here using multiple analysis window lengths allows one to simultaneously detect weak correlations that are persistent in time and stronger correlations that exist transiently. The example shown in Fig. 8 illustrates the failure of the classical coherence estimate for nonstationary data. In this example, the nonstationarity resulted from a shift in the relative phase of the signals causing a reduction in phase coherence.
It is important to keep in mind that significant values only imply that the correlations, as measured by the phase-locking indices, are unlikely to occur by chance. Nonsignificant values do not imply that there are no correlations or no interactions. This is particularly important to consider when the neural interactions under scrutiny might be of very short duration. In these cases, the interactions are likely to be beyond statistical resolution regardless of the method used.
Phase slips in oscillatory data
Within oscillatory segments (i.e., those exhibiting high SNR), the phase representation presents occasional singularities. The singularities that take place in the midst of an oscillatory episode can be interpreted as phase slips, and the amount of phase advance associated with them can be calculated by the method proposed here (Fig. 3). The analysis of phase slips can be useful for understanding the mechanisms whereby oscillatory neuronal populations synchronize or desynchronize. Indeed, phase slips could be the dynamical signature of either a synchronizing or a desynchronizing mechanism. The former could occur for instance if oscillatory cells of disparate natural frequencies become coupled. In this case, phase slips might contribute to keeping the relative phase of the populations within bound.
Apart from their theoretical interest, phase slips are a feature of the dynamics to be incorporated in the generation of surrogate data, particularly in cases where they occur frequently. Indeed, for pairs of independent oscillators, phase slips will tend to accelerate the diffusion of the relative phase, causing a decrease in the phase-locking indices and a subsequent reducing the cutoff values. If slips are not incorporated, the resulting higher cutoff levels are likely to produce false negatives.
To extract slipping events, we took advantage of the fact that phase slips appear as spikes in the instantaneous frequency. However, other situations also give rise to spikes, and a distinction must be made between real phase slips and spikes caused by noise or aperiodic activity. Real phase slips occur during oscillatory epochs and are well defined in the sense that the amount of angular advance is, to some degree, independent of the method used to compute phase. In contrast, spikes that appear during episodes of aperiodic or noisy activity are not flanked by oscillation cycles and occur because of the lack of a well defined phase. For the latter, a phase advance cannot be defined unambiguously. An easy way to distinguish between these two types of spikes is to concurrently use the SNR criterion to determine whether the signal is oscillatory.
Comparison of phase-locking indices
We compared three different phase-locking indices and found, for any surrogate data method, only minor differences in the results. Usually, these differences could be minimized by selecting slightly different cutoff levels for the different methods (say, a 95% cutoff level for the phase coherence could be matched by using a 92% level for the entropy index), yet we did not find any method to be consistently stricter than the others throughout the data sets analyzed. Most often, differences were seen for the shorter analysis windows, which reflect the fact that shorter time windows return less reliable estimates. The differences can also be due to the specific properties of the phase locking utilized. Although phase coherence is a good measure of how clustered the relative phase is around a single value, the entropy measure can detect multimodal clustering (or any deviation from an uniform distribution). The mutual information index has the additional advantage of detecting interdependencies other than a constant phase difference. Most importantly, the mutual information method differs from the other two in that it is computationally more expensive, due to the need for estimating the joint probabilities as a function of time. Using test data generated from simple mathematical models exhibiting synchronous activity may aid in how one would interpret these indices in a biological context.
Concluding remarks
Phase synchrony has been hypothesized to play an important role in many neural systems. The methods presented here can be easily adapted to different types of signals. We have applied it to spike trains and EMG, but it can also be used to study phase correlations in LFP, EEG, and MEG signals. The only essential requirement for this method is that the signals contain a clear oscillatory component. This technique can be adapted to study central pattern generators, some prominent CNS rhythms such as delta sleep (Destexhe et al. 1999), theta rhythms (Harris et al. 2002; O'Keefe and Recce 1993), or the rhythms of epileptic discharges (Neckelmann et al. 1998). However, the method should be used with caution in signals that have broad band spectra or whenever a phase cannot be defined unambiguously.
Statistically significant estimates of episodes of synchrony allow one to draw conclusions about the functional organization of the networks that support them. Applying the method to the particular problem of tremor-generating networks in PD reveals dynamical features that were missed by using standard coherence methods (Fig. 8). In the example of Fig. 8, phase locking could occur episodically between a cell and muscle exhibiting oscillatory activity. The phase entrainment occurred at different lags in both episodes, and this dynamic property was missed by the global coherence estimate (Fig. 8A). Indeed this example shows how nonstationary phase locking can drop the coherence estimate below significance. Using multiple analysis windows enables us to view phase locking as it evolves over time to find intervals where the phase locking is strongest. We are currently using these methods to explore how the phenomenon is distributed within or across structures of the basal ganglia—thalamocortical loop; for example, we are studying the possibility that the parallel circuits exhibiting tremor activity might synchronize episodically, or intermittently. This may help to understand how the organization of the normal motor networks is altered to support tremor generating circuits.
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
GRANTS
This work was supported by National Institute of Neurological Disorders and Stroke Grant NS-39121.
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: K. A. Sigvardt, Center for Neuroscience, University of California, Davis, 1544 Newton Ct., Davis, CA 95616 (E-mail: kasigvardt{at}ucdavis.edu).
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSDISCUSSIONACKNOWLEDGMENTSREFERENCES |
|---|
Boccaletti S, Kurths J, Osipov G, Valladares DL, and Zhou CS. The synchronization of chaotic systems. Phys Reports 366: 1–101, 2002.
Brody CD. Slow covariations in neuronal resting potentials can lead to artefactually fast cross-correlations in their spike trains. J Neurophysiol 80: 3345–3351, 1998.
Brody CD. Correlations without synchrony. Neural Comput 11: 1537–1551, 1999a.[CrossRef][Web of Science][Medline]
Brody CD. Disambiguating different covariation types. Neural Comput 11: 1527–1535, 1999b.[CrossRef][Web of Science][Medline]
Carter GC. Coherence and time delay estimation. Proc IEEE 75: 236–255, 1987.[CrossRef]
Csicsvari J, Jamieson B, Wise KD, and Buzsaki G. Mechanisms of gamma oscillations in the hippocampus of the behaving rat. Neuron 37: 311–322, 2003.[CrossRef][Web of Science][Medline]
Destexhe A, Contreras D, and Steriade M. Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states. J Neurosci 19: 4595–4608, 1999.
Donoghue JP, Sanes JN, Hatsopoulos NG, and Gaal G. Neural discharge and local field potential oscillations in primate motor cortex during voluntary movements. J Neurophysiol 79: 159–173, 1998.
Engel AK, Fries P, and Singer W. Dynamic predictions: oscillations and synchrony in top-down processing. Nat Rev Neurosci 2: 704–716, 2001.[CrossRef][Web of Science][Medline]
Feige B, Aertsen A, and Kristeva-Feige R. Dynamic synchronization between multiple cortical motor areas and muscle activity in phasic voluntary movements. J Neurophysiol 84: 2622–2629, 2000.
Fetz EE. Temporal coding in neural populations? Science 278: 1901–1902, 1997.
Gabor D. Theory of communication. J Inst Elect Eng 93: 429–457, 1946.
Glass L and Mackey M. From Clocks to Chaos. Princeton, NJ: Princeton University Press, 1988.
Gray CM. The temporal correlation hypothesis of visual feature integration: still alive and well. Neuron 24: 31–47, 111–125, 1999.[CrossRef][Web of Science][Medline]
Gray CM, Konig P, Engel AK, and Singer W. Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature 338: 334–337, 1989.[CrossRef][Medline]
Gross J, Tass PA, Salenius S, Hari R, Freund HJ, and Schnitzler A. Cortico-muscular synchronization during isometric muscle contraction in humans as revealed by magnetoencephalography. J Physiol 527: 623–631, 2000.
Harris KD, Henze DA, Hirase H, Leinekugel X, Dragoi G, Czurko A, and Buzsaki G. Spike train dynamics predicts theta-related phase precession in hippocampal pyramidal cells. Nature 417: 738–741, 2002.[CrossRef][Medline]
Hoyer D, Leder U, Hoyer H, Pompe B, Sommer M, and Zwiener U. Mutual information and phase dependencies: measures of reduced nonlinear cardiorespiratory interactions after myocardial infarction. Med Eng Phys 24: 33–43, 2002.[CrossRef][Web of Science][Medline]
Hurtado JM, Gray CM, Tamas LB, and Sigvardt KA. Dynamics of tremor-related oscillations in the human globus pallidus: a single case study. Proc Natl Acad Sci USA 96: 1674–1679, 1999.
Hurtado JM, Lachaux JP, Beckley DJ, Gray CM, and Sigvardt KA. Inter- and intralimb oscillator coupling in parkinsonian tremor. Mov Disord 15: 683–691, 2000.[CrossRef][Web of Science][Medline]
Hutcheon B and Yarom Y. Resonance, oscillation and the intrinsic frequency preferences of neurons. Trends Neurosci 23: 216–222, 2000.[CrossRef][Web of Science][Medline]
Jarvis MR and Mitra PP. Sampling properties of the spectrum and coherency of sequences of action potentials. Neural Comput 13: 717–749, 2001.[CrossRef][Web of Science][Medline]
Kirkwood PA. On the use and interpretation of cross-correlations measurements in the mammalian central nervous system. J Neurosci Methods 1: 107–132, 1979.[CrossRef][Web of Science][Medline]
Konig P. A method for the quantification of synchrony and oscillatory properties of neuronal activity. J Neurosci Methods 54: 31–37, 1994.[CrossRef][Web of Science][Medline]
Lachaux JP, Rodriguez E, Martinerie J, and Varela FJ. Measuring phase synchrony in brain signals. Hum Brain Mapp 8: 194–208, 1999.[CrossRef][Web of Science][Medline]
Le Van Quyen M, Foucher J, Lachaux J, Rodriguez E, Lutz A, Martinerie J, and Varela FJ. Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony. J Neurosci Methods 111: 83–98, 2001.[CrossRef][Web of Science][Medline]
Miller WL and Sigvardt KA. Spectral analysis of oscillatory neural circuits. J Neurosci Methods 80: 113–128, 1998.[CrossRef][Web of Science][Medline]
Mitra PP and Pesaran B. Analysis of dynamic brain imaging data. Biophys J 76: 691–708, 1999.[Web of Science][Medline]
Neckelmann D, Amzica F, and Steriade M. Spike-wave complexes and fast components of cortically generated seizures. III. Synchronizing mechanisms. J Neurophysiol 80: 1480–1494, 1998.
O'Keefe J and Recce ML. Phase relationship between hippocampal place units and the EEG theta rhythm. Hippocampus 3: 317–330, 1993.[CrossRef][Web of Science][Medline]
Palus M and Hoyer D. Detecting nonlinearity and phase synchronization with surrogate data. IEEE Eng Med Biol Mag 17: 40–45, 1998.[CrossRef][Web of Science][Medline]
Pesaran B, Pezaris JS, Sahani M, Mitra PP, and Andersen RA. Temporal structure in neuronal activity during working memory in macaque parietal cortex. Nat Neurosci 5: 805–811, 2002.[CrossRef][Web of Science][Medline]
Pikovsky A, Rosenblum M, and Kurths J. Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge, UK: Cambridge Univ. Press, 2001.
Rodriguez E, George N, Lachaux JP, Martinerie J, Renault B, and Varela FJ. Perception's shadow: long-distance synchronization of human brain activity. Nature 397: 430–433, 1999.[CrossRef][Medline]
Rosenberg JR, Amjad AM, Breeze P, Brillinger DR, and Halliday DM. The Fourier approach to the identification of functional coupling between neuronal spike trains. Prog Biophys Mol Biol 53: 1–31, 1989.[CrossRef][Web of Science][Medline]
Rosenblum MG, Pikovsky AS, and Kurths J. Phase synchronization of chaotic oscillators. Phys Rev Lett 76: 1804–1807, 1996.[CrossRef][Web of Science][Medline]
Rosenblum MG, Pikovsky AS, Kurths J, Schafer C, and Tass PA. Phase synchronization: from theory to data analysis. In: Neuro-Informatics, edited by Moss F and Gielen S. Amsterdam: Elsevier Science, 2001, p. 279–321.
Schafer C, Rosenblum MG, Abel HH, and Kurths J. Synchronization in the human cardiorespiratory system. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 60: 857–870, 1999.[Medline]
Schaltenbrand G and Wahren W. Atlas for Stereotaxy of the Human Brain, With an Accompanying Guide. Stuttgart, Germany: Thieme, 1977.
Schreiber T and Schmitz A. Surrogate time series. Physica D 142: 346–382, 2000.[CrossRef]
Shadlen MN and Movshon JA. Synchrony unbound: a critical evaluation of the temporal binding hypothesis. Neuron 24: 67–77, 111–125, 1999.[CrossRef][Web of Science][Medline]
Singer W and Gray CM. Visual feature integration and the temporal correlation hypothesis. Annu Rev Neurosci 18: 555–586, 1995.[CrossRef][Web of Science][Medline]
Tass P, Rosenblum MG, Weule J, Kurths J, Pikovsky A, Volkmann J, Schnitzler A, and Freund HJ. Detection of n:m phase locking from noisy data: application to magnetoencephalography. Phys Rev Lett 81: 3291–3294, 1998.[CrossRef]
Thomson DJ and Chave AD. Jackknifed error estimates for spectra, coherences and transfer functions. In: Advances in Spetrum Analysis and Array Processing, edited by Haykin S. Englewood Cliffs, NJ: Prentice Hall, 1991, p. 58–113.
Varela F, Lachaux JP, Rodriguez E, and Martinerie J. The brainweb: phase synchronization and large-scale integration. Nat Rev Neurosci 2: 229–239, 2001.[CrossRef][Web of Science][Medline]
This article has been cited by other articles:
|
|
 |
|
  B. Hangya, Z. Borhegyi, N. Szilagyi, T. F. Freund, and V. Varga GABAergic Neurons of the Medial Septum Lead the Hippocampal Network during Theta Activity J. Neurosci., June 24, 2009; 29(25): 8094 - 8102. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 G. A. Jacobson, I. Lev, Y. Yarom, and D. Cohen Invariant phase structure of olivo-cerebellar oscillations and its putative role in temporal pattern generation PNAS, March 3, 2009; 106(9): 3579 - 3584. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 X. Li, D. Cui, P. Jiruska, J. E. Fox, X. Yao, and J. G. R. Jefferys Synchronization Measurement of Multiple Neuronal Populations J Neurophysiol, December 1, 2007; 98(6): 3341 - 3348. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. B. L. Tort, H. G. Rotstein, T. Dugladze, T. Gloveli, and N. J. Kopell On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus PNAS, August 14, 2007; 104(33): 13490 - 13495. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. M. Hurtado, L. L. Rubchinsky, K. A. Sigvardt, V. L. Wheelock, and C. T. E. Pappas Temporal Evolution of Oscillations and Synchrony in GPi/Muscle Pairs in Parkinson's Disease J Neurophysiol, March 1, 2005; 93(3): 1569 - 1584. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
