| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Center for Neuroscience, University of California, Davis; 2Department of Neurology and 3Department of Neurosurgery, University of California Davis Medical Center, Sacramento; and 4Clinical Neuroscience, Kaiser Permanente, Sacramento, California
Submitted 13 August 2004; accepted in final form 2 October 2004
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
60 s), GPi tremor-related units were statistically coherent with restricted regions of the peripheral musculature displaying tremor. The distribution of pooled coherence across all pairs supports a classification of GPi cell/EMG oscillatory pairs into coherent or noncoherent. Analysis using 2-s sliding windows shows that oscillatory activity in both GPi tremor units and muscles occurs intermittently over time. For brain/muscle pairs that are coherent, there is partial overlap in the times of oscillatory activity but, in most cases, no significant correlation between the times of oscillatory subepisodes in the two signals. Phase locking between coherent pairs occurs transiently; however, the phase delay is similar for different phase-locking subepisodes. Noncoherent pairs also show episodes of transient phase locking, but they occurred less frequently, and no preferred phase delay was seen across subepisodes. Tremor oscillations in pallidum and EMGs are punctuated by phase slips, which were classified as synchronizing or desynchronizing depending on their effect on phase locking. In coherent pairs, the incidence of synchronizing slips is higher than desynchronizing slips, whereas no significant difference was seen for noncoherent pairs. The results of this quantitative characterization of parkinsonian tremor provide a foundation for hypotheses about the structure and dynamical functioning of basal ganglia motor control networks involved in tremor generation. |
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
), including several movement disorders (Farmer 2002). Oscillations in the motor cortex and related subcortical structures can propagate to the limb musculature via descending pathways (Baker et al. 2003; Gross et al. 2000; Salenius and Hari 2003) and give rise to tremors (Elble 1996; Elble and Koller 1990; Halliday et al. 2000; McAuley and Marsden 2000). One of the best known is the resting tremor of Parkinson's disease (PD) that affects roughly 70% of PD patients. Parkinsonian tremor occurs in the 2- to 6-Hz frequency band, is predominantly present at rest, and is more prevalent in the distal than proximal segments of the upper or lower limbs (Deuschl et al. 2000).
The pathological origin of PD, the loss of midbrain dopaminergic cells, has been well known for nearly 40 yr (Hornykiewicz 1966), but despite significant advances in our understanding of the anatomy and cellular physiology of the basal ganglia (BG) structures affected by the dopamine (DA) loss, the causal link between the pathology and the symptoms remains obscure. Most recently, hypotheses regarding the pathophysiology of PD have emphasized the dynamical changes that the BG network undergoes as a result of the pathology; in this view, changes in firing patterns, not just of mean firing rates, are considered the hallmark of PD (reviewed in Bevan et al. 2002; Brown 2003; Obeso et al. 2000). This idea has been supported by physiological evidence from PD patients and animal studies, including the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) model of parkinsonism (Nini et al. 1995). It has been shown that following DA loss, cells in globus pallidus internus (GPi) and subthalamic nucleus (STN) enter a regime of rhythmic firing concomitant with excess synchrony (Bergman et al. 1998; Brown et al. 2001; Levy et al. 2002; Obeso et al. 2000). Given the tight interrelationship between the BG and other structures in the motor system, it is not surprising that the pathological regime of PD can affect the entire limb motor network, including cortico-thalamic, cerebellar, and spinal networks (reviewed in Doya 2000; Houk and Wise 1995; Middleton and Strick 2000). Several studies indicate that PD tremor-related activity is not restricted to the core BG/ motor cortex network, but it involves cerebellar and spinal circuitry as well (Parker et al. 1992; Timmermann et al. 2003; Volkmann et al. 1996; Williams et al. 2002).
A prominent feature of PD dynamics is that tremor-related oscillations occurring in different locations of the network can in many cases be uncorrelated, even while their frequencies are nearly uniform. This phenomenon was first studied by Alberts et al. (1969), who compared tremor-frequency intracortical potentials to the EMGs of muscles in different limbs. Several studies of multilimb EMG recordings in PD patients also indicate that tremor in different limbs is largely uncorrelated (Ben-Pazi et al. 2001; Hurtado et al. 2000; O'Suilleabhain and Matsumoto 1998; Raethjen et al. 2000). In addition, dual recordings in GPi tremor-related cells during stereotactic surgery have shown that, although cells might be correlated to restricted portions of the musculature or to each other, uncorrelated oscillations are commonplace as well (Hurtado et al. 1999; Lemstra et al. 1999). These findings are consistent with a topographic organization of the individual structures that comprise the tremor generating network.
However, the results of these studies also raise a number of questions. Clinical observations indicate that tremor in any one muscle can come and go over time; in spectral analysis terms, tremor is largely a nonstationary phenomenon. If nonstationary signals are analyzed over long periods of time, even if they contain short segments of phase locking, the correlation could be below statistical significance. A better assessment of the interactions would be obtained by analyzing these signals over time. Another common observation made during intraoperative mapping of the basal ganglia in functional neurosurgery is that the tremor-related neural activity might wax and wane while the limb tremor persists, and vice versa, persistent tremor-related activity can exist in the apparent absence of limb tremor. This raises a similar and related question regarding the temporal evolution of the correlations between tremor-related neurons within the basal ganglia and between these neurons and the peripheral musculature. Finally, the topographic organization of the different limb representations within the basal ganglia (DeLong et al. 1985) raises the question of how the neurons within and between these representations interact during tremor episodes.
The purpose of this study is to address these questions in detail. We applied time-dependent spectral techniques to tremor-related single unit and EMG recordings in PD patients. Using statistical methods that detect the onset and offset of oscillations and phase locking at high temporal resolution (Hurtado et al. 2004), we characterized the temporal properties of tremor-related oscillatory activity in the basal ganglia-muscle network. The results of this quantitative characterization of tremor dynamics provide a foundation for hypotheses about the structure and dynamical functioning of basal ganglia motor control networks involved in parkinsonian tremor.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
We studied data collected from five cases of surgical treatment of patients with idiopathic Parkinson's disease. This included all patients in our program from 1998 to 2003 who had a microelectrode-guided pallidal procedure (pallidotomy or pallidal deep brain stimulating (DBS) electrode implantation and who fulfilled our inclusion criteria of having recordings of tremor-related activity from the internal segment of the GPi and simultaneously recorded tremulous EMGs (see Table 1). Criteria for identification of neurons as being within GPi included 1) area with noted increased level of tonic activity; 2) area with tremor-related activity; 3) area with passive responses to limb movement; 4) quiet zone below GPi with optic tract below that (as indicated by light responses); 5) lesions placed in area identified in this way result immediately in improvement of one or more motor symptoms (e.g., reduction in tremor and rigidity); 6) identification of correct placement or lesion or DBS in postoperative MRI; and 7) improvement in UPDRS scores 6 mo postoperatively.
|
1 min. Recordings
Data were collected during microelectrode-guided mapping of the posteroventral pallidum prior to lesion or the placement of a DBS electrode. Details of the procedure have been previously described (Hurtado et al. 2004). Briefly, 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 MRIs of the brain in relation to fiducial markings on the stereotactic frame (Radionics, Burlington, MA). EMGs were recorded from one to four of the following arm and leg muscles: abductor pollicis brevis (APB), wrist extensor (WE), wrist flexor (WF), biceps (BI), quadriceps (QUAD), gastrocnemius (GAST), and tibialis anterior (TA). 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 was approved by the Institutional Review Boards of The University of California Davis and Kaiser Permanente Research Foundation.
Microelectrode recordings were obtained with 50-µm, beveled, stainless steel electrode (Frederick Haer and Company, Bowdoinham, ME), with an impedance of 0.5–1.0 M
at 1 KHz. Signals were amplified (10,000 times) and band-pass filtered between 500 Hz and 10 KHz (Bak Electronics, Germantown, MD). EMGs were recorded with Grass disk electrodes (Astro-Med, West Warwick, RI), amplified (2,000 times), band-pass filtered in the range of 30 Hz to 1 KHz, and full wave rectified prior to spectral analysis. This last step was to recover the low-frequency myoelectric signal from the amplitude-modulated activity in the EMG signals. Neuronal and EMG signals were digitized at 10 kHZ using an A/D board (RC Electronics Santa Barbara, CA) with software provided by that company. Neuronal spikes were threshold extracted (>2 SD above baseline) and down sampled at a 1-ms resolution. Single units were further assessed by confirming the presence of a clear refractory period (>1–2 ms) in the autocorrelogram and stability in spike shape and amplitude over the recording period. Multiunit data were discarded from the analysis.
Standard spectral analysis
Standard statistical spectral methods were used to compute overall (i.e., over 60 s) spectral properties of the signals, as described previously (Brillinger 1981; Hurtado et al. 1999; Percival and Walden 1993; Rosenberg et al. 1989). The terms "overall spectra" and "overall coherence" are used to distinguish these computations from the time-dependent coherence measure we apply later. All analysis was performed in MATLAB (MathWorks, Natick, MA). For overall spectral estimation, we used a multiple window method (Carter 1987; Percival and Walden 1993) as follows: epochs of paired neural-EMG recordings were subdivided into Nwin 1-s segments with no overlap and multiplied with a Hanning tapering window to reduce spectral leakage. The discrete Fourier transform (DFT) of each segment di,k (channel i, segment k) was computed using the FFT algorithm from the MATLAB library. The frequency bin for the DFT and all overall spectral estimates are 1 Hz. To obtain the power spectrum, a direct spectral estimate (i.e., periodogram) from each 1-s segmentwas computed by taking the magnitude square of its DFT ateach frequency, 
where 
is the frequency, andthe spectrum is estimated as the average of all direct estimates, 
To extract significant peaks, we consider as a null hypothesis a flat spectrum with the same total power as the observed spectrum as follows. For a stationary process, the spectrum follows a 
2 distribution (Brillinger 1981; Percival and Walden 1993)
 | (1) |
null is the empirical (i.e., sample) spectrum in the null condition, and E[Pnull] is the expected value. For Nwin-independent windows, the degrees of freedom are v = 2N. Under the null hypothesis of a flat spectrum, E[Pnull] should be constant across frequencies. For evaluation purposes, the constant value was taken as the average of the empirical spectra in 
. This assures that the null distribution has equal power as the empirical spectrum. In our case, the average was taken in the frequency range between 1 and 35 Hz. Departures from the null condition (i.e., peaks) are obtained at 
= 0.002 significance level, and significant peaks in the range of tremor (2–6 Hz) were considered as evidence of tremor-related activity.
Overall coherence is computed from the formula
 | (2) |
ij() is the estimated cross spectrum between channels 
, where the bar denotes complex conjugate, and Nwin is the number of nonoverlapping analysis windows.
To distinguish between noncoherent and coherent pairs, we test the null hypothesis of independent activity. We define as a cut-off level 
as the value of coherence such that the probability P( > |
= 0) = . This is the maximum value consistent with the zero coherence hypothesis at the significance level . A cut-off value as a function of and N is given by
 | (3) |
). The validity of this binary classification (coherent/noncoherent) was tested by building a histogram of coherence values across all recorded pairs and checking the modes. The method is described below. An example of a digitized EMG and threshold-extracted spike train is shown in Fig. 1A. The overall power spectrum of each of these signals and their coherence spectrum is shown in Fig. 1B. Recording pairs having a significant peak in their coherence spectrum are expected to maintain a relatively constant phase lag at that frequency throughout the recording period. We display the distribution of phase differences over the sample windows as an angular histogram (Fig. 1B). Clustering around an angular bin is expected when coherence is high.
|
Building a histogram of combined or "pooled" coherence values across patients, neural/EMG pairs, and recording episodes is complicated by two problems. First, coherence estimators obtained with different numbers of independent windows (N) are not comparable since the bias and variance of the estimator depends crucially on N (Amjad et al. 1997; Carter 1987). Second, our number of recorded pairs was too limited to generate a meaningful histogram. We adopted a bootstrapping strategy to address these problems. In this scheme, the following sequence was repeated n = 105 times: first, a recording session was selected at random from the database, with replacement; the probability of choosing any particular session was made proportional to that session duration. From the chosen session a neural/EMG pair is obtained; if more than one EMG channel was available, one was selected randomly. Second, the paired neural/EMG data were subdivided into nonoverlapping 1-s segments, as described in the previous section. Third, from these (paired) segments, a group of 20 is randomly drawn, and their coherence is computed. The total of 105 coherence values obtained was pooled to form a histogram; pooling was legitimate since each coherence estimate was computed from the same number of segments (n = 20).
The bias and variance of the coherence estimator depends also on the "underlying" coherence (Carter 1987). To equalize variances across coherence values, we used the Fisher transform of coherence, 
, which returns values close to a normal distribution with variance that is not dependent on the underlying coherence (Amjad et al. 1997).
Time-frequency analysis, SNR, and tremor-on-off index
We implemented a method to extract the intervals of significant oscillatory activity in the EMG and neuronal data. This step is important for the phase reconstruction procedure that follows, since a phase can be meaningfully defined only for signals that are oscillatory. First, a time-frequency spectrum, (,t) (, frequency; t, time) was computed for each data series. For each point in time, we applied the same calculation scheme as for the overall spectral analysis, with the sole difference being that here we used overlapping windows (90%) to obtain a higher temporal resolution. Data were windowed into 1-s sliding windows (100-ms offset between subsequent windows), a Hanning taper was applied to each window, and the direct spectra was computed from its FFT, as described before. The spectra of 10 adjacent windows were averaged to yield a spectral estimate. One such estimate was obtained every 100 ms, spanning 1.9 s of data (i.e., 10 1-s windows 100-ms offset) and with a frequency resolution of 1 Hz.
Episodes of tremor activity are defined as intervals where the spectral peak in the tremor range (2–6 Hz) is substantially higher than the baseline (noise) level. A cut-off level was obtained bytaking a flat power spectrum as a null hypothesis and considering departures at a significance level = 0.002 as significant peaks. The null distribution of the averaged spectrum, null, is 2, and we approximate E[Pnull] by the average of the power in the frequency range 1–35 Hz, . For the degrees of freedom, we must consider that, in this case, there is substantial overlap between adjacent windows, and the assumption of independence is thus violated. We therefore estimated the equivalent degrees of freedom for nonindependent windows from a procedure described in detail in Thomson and Chave (1991). For an overlap level of 90%, usingNwin =[r] 10 and considering Hanning tapers, the degrees of freedom parameter is v = 4.87. At a significance level = 0.002, the corresponding cut-off value is 3.78 
. Peaks that surpass the spectral average baseline 3.78 times were thus considered significant (P < 0.002). Since we are extracting phases in the tremor range, we consider peaks in the 2- to 6-Hz range.
To implement this criterion over time, we defined a signal to noise ratio (SNR) as a function of time, SNR(t), as the ratio of maximum power in the tremor frequency band over the average signal power in the broader band 1–35 Hz
 | (4) |
a, b are the limits of the tremor band (2–6 Hz), max = 35 Hz, and 0 = 1 Hz. From the previous considerations, we used a threshold value of SNR(t) > 3.78 to extract the segments where tremor oscillations are significant. We confirmed the validity of the threshold by visually inspecting tremor-related neural activity and tremor EMG series and checking that this threshold index discriminates well between oscillatory and nonoscillatory intervals (Fig. 1C). Intervals above and below this threshold are referred to as "tremor-on" and "tremor-off" intervals, respectively. It is important to note that this separation between nonoscillatory and oscillatory episodes is necessary for a valid phase reconstruction, since phase variables defined over nonoscillatory data are misleading.
Tremor-on and tremor-off intervals were represented by a binary variable as a function of time 
(t), as
 | (5) |
= 0 for tremor-off intervals and = 1 for tremor-on intervals). This index was computed for neuronal activity and EMG recordings. Since one SNR point, spanning 1.9 s of data, was obtained every 100 ms, gaps in tremor (i.e., periods of = 0) that lasted < 1 s and were in the midst of tremor-on periods ( = 1), were considered nonsignificant and set to the value = 1. We did this correction because a drop in the SNR value for such short times does not reflect a lasting interruption in tremor (the total time extent for computing a SNR point was 1.9 s, i.e., 10 1-s windows with 90% overlap). These short interruptions were caused by the appearance of a series of phase slips in quick succession and are treated separately in phase slips and estimation of associated phase advance. This transformation from continuous SNR to binary indexes provides a way to extract episodes of significant tremor, where a phase can be defined with relative certainty, and it simplifies the analysis of temporal incidence of tremor.
Tremor-on-off correlation for paired recordings
To determine whether there are temporal correlations between the tremor-on-off intervals for different neural and EMG channel pairs, we computed a (0 lag) correlation coefficient Ri,j between the of channel pairs (i,j), over each recording episode
 | (6) |
close to 1 means that channels i,j tend to be in the tremor-on or -off states concurrently, whereas a value close to 0 indicates that the tremor-on and -off states in both channels are independent; a value near –1 means that one channel is likely to be tremor-off when the other is tremor-on. The statistical significance of Ri,j was assessed by estimating the distribution of correlation coefficients under the null hypothesis of 
and obtaining cut-off values from this distribution. To estimate the distribution, independent pairs of i(t) and j(t) series were generated using surrogate data. A set of surrogate series is(t) can be constructed by time-shifting either of the original series by a random step, while keeping circular boundary conditions. The Ri,j,s (superscript s for surrogate) was computed for 500 surrogate data pairs, and the 95th and 5th percentiles of the distribution were used as the cut-off to determine significant correlation (or anticorrelation). In cases where either of the i,j(t) was 1 or 0 for all t, the variance is zero and Ri,j has no meaning. A caveat of this analysis is that there may be a loss of information when transforming the SNR measure to a binary index, and therefore significant correlations in SNR could be masked. We therefore computed the correlation coefficients of the original SNR, and the answers were identical to those obtained by using the tremor index.
Phase reconstruction procedure
The goal of the phase reconstruction procedure is to obtain an angular phase variable as a function of time from the oscillatory time series. The phase variable indicates the amount of completion of an oscillation cycle at any given point in time and is used to quantify phase locking between oscillatory signals. Importantly, a phase reconstruction method yields meaningful results only if the signal is oscillatory (Hurtado et al. 2004; Pikovsky et al. 2001; Rosenblum et al. 2001). We therefore restrict the analysis of phase to tremor-on intervals.
The steps in the phase reconstruction methods have been described in detail elsewhere (Hurtado et al. 2004). In summary, after spike extraction in the neuronal signal and rectification in EMGs, signals were band-pass filtered in the range of tremor activity. We used a digital FIR filter, 2- to 6-Hz passband, stop 0–1 (–60 dB) and 7 Hz Nyquist, (–80 dB), sampled at 1 kHz. The parameters yielded an almost flat response in the passband, and rejection of out of band activity with a sharp rolloff. A filter of relatively high order (2,530) was necessary to obtain good out of band rejection together with a flat response in the tremor range and little distortion in the time location of peaks. We compared the original series to the filtered ones to check that filtering did not alter peak locations and that it was sensitive enough to follow the shifts in phase in the original data. The process of filtering entails a loss of temporal resolution, since it is equivalent to convolving the data with a moving window. High filter orders such those used here imply longer windows. We measured the time as localization of our filter by taking the square of its impulse response as a distribution (Percival and Walden 1993) and computing its SD. In our case, this was 188 ms, meaning that the lost of temporal precision due to filtering happened within a very narrow window of time, even though the total filter length was about 5 s (Fig. 1D).
Phase is obtained from the filtered series by using the Gabor representation, 
(t), that projects the oscillatory series on to the complex plane and results in a rotational trajectory around the complex origin (Gabor 1946; Rosenblum et al. 2001). Its real part is the signal itself Re[(t)] = x(t), and the imaginary part is obtained from the Hilbert transform of the signal Im[(t)] = 
(t) = H[x(t)]. The complex trajectory (t) is projected onto the unit circle (Fig. 1E): 
, and the phase of the signal is obtained as the angular variable 
. As a result, we obtained a discrete series of 
(tk), k = 1, 2...N, where 
t = tk+1 – tk is the sampling period (1 ms), and N is the number of samples. The "raw" instantaneous frequency (t) is defined here as the first-order differential of 
. Ideally, varies slowly over time (Boashash 1992); the problem of discontinuities is treated in Phase slips and estimation of associated phase advance.
Phase slips and estimation of associated phase advance
Phase slips are discontinuities in the phase evolution of an oscillator, leading to abrupt phase advances. To detect slips, we took advantage of the fact that they would appear as "spikes" in the derivative of phase, which is the instantaneous frequency, (t) (Fig. 2A, arrows). We therefore referred to these events as -spikes. -spikes were extracted from (t) by setting as threshold the frequency limits of the band-pass filter that was used in the phase construction procedure (Fig. 2A, dashed lines). The rationale for this choice is that the instantaneous frequency is expected to remain within the frequency band of the signal (Boashash 1992); large departures can be considered singularities. A flat band-pass response and good out of band rejection of the filter were important in the procedure of phase slip extraction because they guaranteed that the only out of band activity visible corresponded to singularities (i.e., phase slips) and minimized false positives in slip detection due to poor filtering. Importantly, only slips occurring in the midst of an oscillation episode, (i.e., with high SNR) were considered as such.
|
(Fig. 2B). The smoothed slips were first "erased" from the raw instantaneous frequency approximation, and the gaps (threshold crossing ± 100 ms) were filled by cubic interpolation from the neighboring endpoints. was obtained by filtering the resulting series (moving average filter, 750 ms width). The forecasted phase series, 
(t) was computed as the time integral of . The additional phase advance due to the slip is the difference between the forecasted phase (with no slip) and the observed phases, 
.
Importantly, this estimate has a margin of error due to the smoothing procedure used to obtain the forecasted phase advance, . The error arises because is the integral of the smoothed instantaneous frequency , whereas equals the integral of the raw instantaneous frequency (see Fig. 2B). Even if no phase slip takes place, the smoothed and raw series would differ, and 
would drift away from the zero expected value. The error should be accounted for in the estimate of for phase slips. We studied the distribution of that error by computing for randomly selected data segments where no slips occurred. From this distribution (data not shown), we found that the angular departure was || < 45° in 95% of the cases. Consequently, slips less than ±45° in magnitude were considered nonsignificant, and a conservative error bound of ±45° was added to all slip angle estimates.
To determine whether there is a preferred angle of phase slipping for a given channel, we studied the distribution of slip advance angles for each recording session and tested the null hypothesis of a uniform circular distribution. An eight-bin circular histogram (45° bin size) of phase advance values was built, and an entropy index was calculated from this histogram. To obtain the null distribution of the entropy index, we numerically computed the index for uniformly distributed angles in the circle. An identical number of events was generated, and the index was computed a 100 different times. From this null distribution, a cut-off value was obtained at the 5th lowest percentile. Entropy values less than the cut-off were indicative of significant clustering.
Phase correlation statistics
To quantify time-locking between pairs over time, we used a time-dependent phase coherence index (details in Hurtado et al. 2004). The phase coherence is always 
1, taking a value of 1 only when the relative phase 
remained constant throughout the observation period. A similar measure has been used by others (Lachaux et al. 1999; Rosenblum et al. 2001). For each data pair, the relative phase series is obtained by taking the difference of angles between the individual phase series 
where the tj is the sampling point (see Fig. 1E).
The time-dependent phase coherence is defined as 
, where N indicates the number of consecutive data samples to be considered in the coherence computation. A small value of N provides an index with a higher temporal resolution, but lower statistical power. In this study we used n = 1,500, which at the sampling rate of 1 kHz corresponds roughly to six cycles of tremor oscillation. Values of phase coherence were only considered when both channels had significant oscillations (i.e., SNR above threshold).
We used a statistical method to test for phase locking against the null hypothesis of independent processes by studying the distribution of for independent processes with similar temporal characteristics ind (see Hurtado et al. 2004). Briefly, surrogate versions of the instantaneous frequency, 
, are created from . To do this, we created surrogates having the same power spectrum as the series, randomizing the phase spectrum while preserving the amplitude spectrum. We erased the -spikes caused by phase slips prior to the computation of surrogates because -spikes contribute to power at all frequencies. From the surrogate , we obtained a series of surrogate 
that are entered into the computation of the null distribution ind. From ind, we obtained a 95% cut-off level ind(.95), to detect times of significant synchronization. In addition, we used the 50% cut-off level, ind(.5), as a reference to obtain the duration of intervals as described below. From the previous analysis step, we could obtain periods of phase locking at a 0.05 significance level. However, even for independent processes, some periods with > ind(.95) above the 95% cut-off ( > ind(.95)) will randomly occur (on average, 5% of the time). To determine whether events above cut-off actually reflect a significant transient synchronization, we examined their duration. For a set of 200 surrogate independent pairs (generated as specified above), we determined the distribution of duration for intervals above the 50% cut-off level ind > ind(.5). A recorded pair was considered to show transient coherence if 1) it contained at least one interval > ind(.5) with duration in the upper 95% and 2) if the within the interval > ind(.5) in question crossed the ind(.95) line at least once.
Synchronizing and desynchronizing phase slips
We describe here a statistical method to determine whether a phase slip in a channel promotes or destroys the synchrony between it and a second (reference) oscillatory signal (Fig. 2C). For this analysis, only phase slips that occur during an interval of synchronization between the two channels are selected. A time segment around the -spike associated with the slip is selected (tsl ± 850 ms, where tsl is the time of threshold crossing) and the phase coherence around the slip occurrence, sl, is computed over that interval. For the computation of sl, the interval tsl ± 100 ms is eliminated from the calculation because the particular shape of an -spike is not a robust feature of the dynamics but depends on the details of numerical methods (filter settings, etc.). A random phase advance chosen from an uniform distribution in the circle is introduced into the channel with the phase slip at t = tsl and the value of with the random slip, R__sl, is computed. This procedure is repeated 100 times, and the histogram of R__sl is computed. A slip is considered "synchronizing" if its sl is in the highest 25th percentile of the R__sl distribution, "desynchronizing" if it is below the lowest 25th percentile, and "neutral" if it falls within the median 50th percentile (Fig. 2D). This procedure was carried out for all single unit recordings, taking each muscle as a reference (neuronal slips), and also for each muscle, taking the neuronal signal as a reference (EMG slips).
A bootstrapping test was performed to determine whether synchronizing slips are significantly more probable in pairs that are coherent overall versus pairs that are noncoherent. Slips from all recordings in our database with overall coherence were listed, and groups of 30 were successively drawn randomly and with replacement. For each group, the number of synchronizing, desynchronizing, and neutral slips was counted. This procedure was repeated 5,000 times, drawing a new set of 30 slips each time. The frequencies of each event type (synchronizing, desynchronizing, and neutral) were averaged over all repetitions to yield an estimated probability distribution for each case. The same procedure was done for slips occurring in noncoherent pairs, and also for slips in EMGs.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
The second set of measures characterizes the temporal evolution of oscillatory activity and the synchrony between oscillations over time, with a resolution of about six tremor cycles (i.e., 1.5 s). From them, we obtain 1) a time profile of the intervals during which tremor-range oscillatory activity is present (tremor-on) and absent (tremor-off) in each of the channels, 2) the time location of phase slips in each of the oscillatory channels during the tremor-on states, and 3) a picture of the time evolution of phase locking between pairs of oscillatory channels. These high resolution measures resolve the problems of nonstationarity faced in the standard spectral measures used to compute overall coherence, that span longer analysis periods.
Distribution of overall tremor coherence reveals two classes of GPi/EMG interactions
Power spectra were computed for GPi single units (n = 9) and EMG recordings (n = 27). The duration of data epochs varied from 30 to 90 s. All single units studied exhibited a statistically significant peak in the tremor range (2–6 Hz), and 22 of the EMG episodes recorded had significant tremor. The resulting 22 neural/EMG pairs with prominent tremor-range activity in both channels amount to a total of nearly 20 min of paired tremor data from all patients, pairs, and recording epochs.
Figure 3 shows the pooled coherence histogram from these 22 pairs, constructed with the bootstrapping method described previously. A total of 105 coherence values was computed from data segments drawn randomly from the database, as described in METHODS. The histogram exhibits a prominent mode at low coherence (arrow) and a wider peak at higher coherence values that contains two closely spaced secondary peaks. The Fisher transform (inset), normalizes the coherence measure and reveals the modes more clearly.
|
The presence of a distinct population of independent oscillators (i.e., 0 coherence) led us to sort our data into coherent and noncoherent pairs; the main motivation for this was to compare the temporal dynamics of one class versus the other. To this end, we used the following classification scheme: a recording was classified as noncoherent if the null hypothesis of independent oscillators could not be rejected at the 0.01 significance level. Pairs that did not conform to the null hypothesis were classified as coherent. We used the full length of the each data epoch to compute coherence for this classification. In addition, the same bootstrapping method was applied, separately, to both coherent and noncoherent subclasses to obtain their individual distributions of coherence values. The black full and dashed lines in the histogram plot show the coherence distributions obtained for these two subsets. The coincidence between the peaks on the latter and those of the full distribution indicates that our classification scheme sensibly splits the population into natural modal components. The full distribution is, of course, the sum of both groups. Note that the distribution of coherence for cell/muscle pairs classified as noncoherent (full line) is largely coincident with the theoretical distribution for independent (i.e., 0 coherence) white noise processes (white line), as expected.
The incidence of coherent and noncoherent pairs is shown at the bottom of Table 2. Of the 22 GPi/EMG pairs that had concurrent oscillations in the tremor range, 17 pairs (77%) showed significant coherence; the remaining 5 (23%) pairs were noncoherent. That many of the pairs were coherent overall is expected, given that both the unit and the muscle in each pair had tremor-related activity. The finding of concurrent, yet noncoherent, oscillatory activity at the same frequency is, however, not intuitive. One possibility is that a given tremor GPi cell is functionally linked in a preferential way to a subset of muscles, but is independent from the rest. This interpretation was supported by recordings where GPi/ tremor correlations were clearly specific to some tremulous muscles but not others. In four of the nine single units studied, we observed overall coherence between an oscillatory GPi unit and some of the recorded muscles but not others, even as these others had significant tremor. In one of our patients (patient B), we recorded two GPi units and four arm muscles. One of the units was coherent with all four muscles, but the other unit was coherent only with the three more distal muscles (APB, WE, WF), but not with the proximal muscle (BI).
The idea that tremor/GPi correlations follow a topographic order was further supported by the results of the analysis of recordings from patients with both upper and lower limb tremor. One of the patients studied here (Table 2, C) exhibited prominent tremor in both upper and lower extremities contralateral to GPi. Overall coherence between a GPi single unit and EMGs was significant only for the two leg muscles, but not for the arm muscles (Fig. 4). This limb selectivity occurs despite the fact that the GPi single unit and all four muscles (2 arm and 2 leg) recorded had prominent oscillatory activity at 4 Hz, as visible in the power spectra. This observation complements our previous finding (Hurtado et al. 1999) that GPi oscillators during tremor are clustered according to limb preference and is in agreement with previous observations that GPi neurons have kinesthetic responses that are limb specific (DeLong et al. 1985).
|
Recorded sessions with significant tremor activity (i.e., a significant peak in the power spectrum) were studied to determine how tremor activity in GPi and EMG evolves over time. The tremor-on-off index, which provides information on the onset and offset of significant tremor oscillations, was used for this purpose. A basic observation is that an episode of tremor-related activity (as well as muscle tremor) characterized by a significant peak in the power spectrum, consists of several tremor-on subepisodes interleaved with tremor-off lapses (Fig. 5, A and B ).
|
; Percival and Walden 1993), shows a large drop in EMG tremor power in 1 compared with 2, with relatively unchanged power between the two episodes in GPi (Fig. 5E). This example shows the situation in a coherent GPi/muscle pair where tremor-related activity persists in GPi as tremor disappears and reappears in the EMG, although the converse situation, where tremor activity persists in the EMG but is intermittent in GPi, was also common.
How is it possible for an intermittent oscillatory pair, where there is only partial overlap in subepisodes of tremor activity, to have statistically significant overall coherence? Inspection of the phase lag histogram for the entire time period (Fig. 5C; see METHODS) indicates that there is significant clustering of the phase delay at 90°; whenever there is concurrent tremor activity, the phase difference returns to a near constant angle. When analyzed over the entire period, this behavior gives rise to high coherence.
We also performed the SNR correlation analysis for muscle/muscle pairs. Of 17 muscle pairs that were coherent at the tremor frequency, the occurrence of tremor was significantly correlated in 7 (41%) and was not in 10 (59%). This analysis of muscle/muscle pairs reveals similar independent tremor-on behavior that was found for the GPi/muscle pairs, even for muscle pairs within the same limb or limb segment.
Phase-locking can be transient in both coherent and noncoherent GPi-muscle pairs
The phase reconstruction outlined in METHODS allows us to track the evolution of phase locking between oscillatory neuronal and muscle activity over time. Since phase reconstruction is meaningful only when oscillatory activity is present, this analysis was restricted to intervals where both channels were in the tremor-on state. In Fig. 6, we show the time dependence of phase coherence for two channels of the same recording session as presented in Fig. 4. Several features are revealed in these plots. First, although GPi and EMG signals are clearly oscillatory for most of the recording time, statistically significant phase locking is episodic and occupies only a fraction of the time that the two signals are concurrently oscillatory. Figure 6A shows the time evolution of phase coherence for a GPi/gastrocnemius pair with overall coherence. There are two periods of concurrent oscillatory activity (a, b); in both of them, significant phase locking is achieved for only a fraction of the time. In period a, significant phase locking begins about 4 s after the onset of concurrent oscillatory activity and ends when tremor stops in the leg muscle. In period b, phase locking also begins after the onset of concurrent tremor. Figure 6B shows the time evolution of phase coherence for a noncoherent pair of the same recording of GPi as Fig. 6A but with respect to an arm muscle (wrist extensor). This pair is statistically noncoherent overall even though tremor activity in the unit and the muscle is concurrent throughout the epoch. Despite being noncoherent, we do, in fact, find episodes of significant transient phase locking, but these episodes are of shorter duration than in the coherent pair. In summary, transient coherence is common in both pairs that are coherent overall and noncoherent pairs. In the analysis of all pairs (see Table 2, transient coherence column), we found that 20 of 22 pairs that had tremor in both channels showed at least one transient coherent episode of significant duration (see METHODS for significance criterion) and 3 of the 5 noncoherent pairs showed significant transient episodes.
|
Close examination of the time course of oscillatory activity reveals the presence of phase slips; these events are characterized by a sudden shift in phase, usually coincident with a decrease in oscillation amplitude lasting less than one tremor cycle (see Fig. 7). Phase slips occur relatively frequently in brain (Fig. 7, A–C) and in muscle (Fig. 7D). We determined the rate of occurrence of these events in GPi and in muscle during each tremor episode (Table 3) and estimated the amount of phase advance associated with each event (see METHODS). Figure 7B shows an example of a phase slip of 190°. If this cell is synchronized to others in the basal ganglia-thalamocortical network, this phase slip would cause a dramatic change in its phase delay with respect to the rest of the network. However, as we see in Fig. 7C, a series of closely spaced phase slips can sum so that the combined phase advance is zero.
|
|
Figure 6, A and B, shows examples of synchronizing and desynchronizing phase slips (1 of each is presented at higher temporal resolution in Fig. 8). Examination of the temporal evolution of coherence in Fig. 6A shows that there is a high incidence of synchronizing phase slips in this coherent pair. In contrast, most of the phase slips in GPi in the noncoherent pair in Fig. 6B occur at points in time where there is no phase locking and therefore have no effect in the already random phase delay. During the periods of transient phase locking, we see one synchronizing and one desynchronizing phase slip.
|
|
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
In this study, we combined classical spectral analysis with a statistical time-dependent method to scrutinize the dynamics of parkinsonian tremor, one of the landmark symptoms of basal ganglia pathophysiology in PD. Our data base included only 9 single units and 27 single unit/EMG pairs. This data
, to normalize the data. Black lines are the bootstrapped coherence obtained from pairs classified as coherent overall (dashed) and noncoherent (solid). White line is the coherence distribution obtained by applying an identical procedure to a numerical model of 2 independent white noise processes; number of bootstrapping samples in the model is chosen to be equal to the number of noncoherent draws (black line).
