| 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
