## Abstract

Synchronous oscillatory dynamics in the beta frequency band is a characteristic feature of neuronal activity of basal ganglia in Parkinson's disease and is hypothesized to be related to the disease's hypokinetic symptoms. This study explores the temporal structure of this synchronization during episodes of oscillatory beta-band activity. Phase synchronization (phase locking) between extracellular units and local field potentials (LFPs) from the subthalamic nucleus (STN) of parkinsonian patients is analyzed here at a high temporal resolution. We use methods of nonlinear dynamics theory to construct first-return maps for the phases of oscillations and quantify their dynamics. Synchronous episodes are interrupted by less synchronous episodes in an irregular yet structured manner. We estimate probabilities for different kinds of these “desynchronization events.” There is a dominance of relatively frequent yet very brief desynchronization events with the most likely desynchronization lasting for about one cycle of oscillations. The chances of longer desynchronization events decrease with their duration. The observed synchronization may primarily reflect the relationship between synaptic input to STN and somatic/axonal output from STN at rest. The intermittent, transient character of synchrony even on very short time scales may reflect the possibility for the basal ganglia to carry out some informational function even in the parkinsonian state. The dominance of short desynchronization events suggests that even though the synchronization in parkinsonian basal ganglia is fragile enough to be frequently destabilized, it has the ability to reestablish itself very quickly.

## INTRODUCTION

The functional significance of oscillations and synchronization in the brain has been extensively studied, for example in relation to perception and cognition (Buzsaki and Draguhn 2004; Engel et al. 2001) and to movement (Baker et al. 1999; Murthy and Fetz 1996; Sanes and Donoghue 1993). Excessively strong, weak or otherwise improperly organized oscillatory activity may contribute to the generation of symptoms of different diseases (Schnitzler and Gross 2005; Uhlhaas and Singer 2006). Oscillations in basal ganglia (BG), in particular at rest and in relation to motor control, are observed in different mammals (rodents, monkeys, humans), in different behavioral conditions, and in different dopaminergic states, e.g., Parkinson's disease (PD) versus normal (Boraud et al. 2005; Gatev et al. 2006). Overall the low-dopamine state as is seen in PD is marked by an increase of oscillatory and synchronous activity (Bergman et al. 1998; Dejean et al. 2008; Priori et al. 2004; Soares et al. 2004).

Movement attenuates the power of oscillations in and synchronization between subthalamic nucleus (STN)/pallidal local field potentials (LFPs) and cortical electroencephalograms (EEG) (Cassidy et al. 2002) and attenuates synchronization between STN neurons (Foffani et al. 2005; Levy et al. 2002). The strength of beta oscillations in STN LFP and in single units is inversely correlated with motor performance in PD (Amirnovin et al. 2004; Kühn et al. 2004). Beta-band synchronization of the EEG of motor areas is correlated with the severity of motor symptoms and decreases as the symptoms were alleviated by different therapies (Silberstein et al. 2005). These and other studies give rise to a “pro-kinetic gamma and anti-kinetic beta” paradigm and underscore the importance of beta-band synchronization for the hypokinetic symptoms of PD (Brown 2007; Hammond et al. 2007; Hutchison et al. 2004).

Correlations of this oscillatory activity between different locations in BG-thalamocortical circuits depend on the brain state (Magill et al. 2006). The spatial organization of synchrony is complex (Goldberg et al. 2004; Lalo et al. 2008). The cellular and network mechanisms of this oscillatory activity are not known; neither are the properties of the temporal dynamics of this synchronization understood. However, understanding the dynamical nature of this synchronization is essential for the understanding of its function as well as for determining potentially efficient therapeutic ways to suppress this synchrony in PD.

This study explores this dynamical nature of synchronous beta oscillations in PD. We hypothesize that this synchrony is strongly variable even on short time scales and investigate how and in what manner synchronous activity is interrupted by asynchronous events. The fragile temporal structure of this synchrony necessitates the use of time-sensitive data analysis methods, which are gaining a wider usage in neuroscience (Le Van Quyen and Bragin 2007). We simultaneously record spikes and LFP from STN of patients with PD and analyze the phase locking between these signals as it develops in time to explore the variability of beta-band oscillations in BG. Phase locking with some (not necessarily zero) phase difference indicates the presence of temporal coordination between signals, and this is reflective of synchrony in its broad sense (Pikovsky et al. 2001). Recorded LFP may primarily reflect synaptic/dendritic activity, while extracellular spikes reflect somatic/axonal activity, thus the phase-locking we study reflects the synchronization between these two modalities.

## METHODS

### Patients and surgery

Eight patients with Parkinson's disease [5 male; age, 63 ± 7 (SD) yr; disease duration, 11 ± 5 yr], who underwent microelectrode-guided implantation of deep brain electrodes in the STN in Indiana University Hospital were included in the study. These are all patients for whom good recordings in STN were available. All patients exhibited hypokinetic symptoms and exhibited no or only weak rest tremor. The decision to perform the surgery was not influenced by subsequent inclusion of the data in the present study. The surgical procedure was carried out using intravenous sedation with dexmedetomidine and local anesthesia and followed standard stereotactic surgical protocols as described in the literature (Machado et al. 2006). Patients were placed in a Leksell stereotactic frame and a contrasted volume-acquisition MRI scan was performed and transferred to a Stealth Station intraoperative navigation computer (Medtronic Navigation, Louisville, CO). Target selection was performed indirectly with respect to anterior commissure (AC) and posterior commissure (PC). Preliminary targets in left and right STN, respectively, were selected with coordinates: *x* = 12.0 mm left (right) of midline; *y* = 2.0 mm behind midcommisural point; *z* = 4.0 mm below AC-PC plane. Coordinates could be modified slightly to accommodate individual variations in MRI-documented anatomy. An FHC microdrive device (FHC, Bowdoin, ME) was used to advance a microelectrode toward and past the previously selected preliminary target. The microelectrode was considered to be within STN when dense, somewhat irregular, high-amplitude action potentials were suddenly recorded after a period of relative silence as the microelectrode passed through the fields of Forel and zona incerta below the thalamus. Finally, the DBS electrode was implanted and its location was confirmed by postoperative MRI. All patients were at least on levodopa and carbidopa containing medications. At the time of surgery, patients had been off antiparkinsonian medication for ≥12 h. The protocol was approved by Indiana University IRB.

### Data recording and processing

Recordings were obtained with 80% platinum/20% iridium glass-insulated microelectrodes (FHC), with an impedance, measured in the brain at 1 kHz, being in the range of 0.5–1.0 MΩ. The recordings were made with Guideline System 3000 (Axon Instruments, Foster City, CA). The recording system amplified the signal (5,000 times) and filtered it into two frequency bands: 300 Hz to 5 kHz and 0–200 Hz to obtain spiking neuronal units and LFP, respectively. Spiking extracellular activity signal and LFP were digitized at 20 kHz rate and saved for off-line analysis.

### Data analysis

Neuronal spikes were threshold extracted (>3 SD above baseline), and single units were isolated with SciWorks/Experimenter software (DataWave Technologies, Berthoud, CO). Both spikes and LFP signals were further downsampled at a 1 kHz rate. Time-series analysis was done with the programs custom written in MATLAB (Mathworks, Natick, MA).

The spectral band of 10–30 Hz was band-pass filtered in both extracellular spiking activity [presented as a binary 0 or 1 (no spike or spike) signal] and in LFP to study oscillatory activity in the beta band. Although the spiking signal was high-pass filtered earlier to detect the times of spikes, filtering to the 10–30 Hz allows us to detect the modulation of the firing rate (bursting in the beta-band). We used an order 1117 Kaiser windowed digital FIR filter sampled at 1 kHz; pass-band 10–30 Hz with pass-band ripple 5%, stop-band 0–8 Hz (−40 dB) and 32 Hz-Nyquist (−40 dB). Although this kind of a filter implies a convolution with a relatively long time-window, there was not much of a temporal information loss because the “temporal localization” of the filter (measured as the SD of the square of its impulse response) was 36 ms. Zero-phase filtering (filtering in both directions) was implemented to avoid phase distortions. The exact spectral boundaries of beta-band activity in the basal ganglia are not known (and may not exist), so we took a relatively broad spectral window, which should accommodate the activity of interest. To detect the intervals of the oscillatory beta-band activity, the following signal to noise ratio (SNR) criterion was implemented. SNR was defined as _{ωa≤ω≤ωb}{P(ω, t)} is the maximum value of the spectral power P(ω, t) in the spectral band of interest [ω_{a}, ω_{b}] and avgω_{0}≤ω≤ω_{max}{P(ω, t)} is the average value in the broader band [ω_{0}, ω_{max}] to include most of the signal. P(ω,t) is obtained by computing the discrete Fourier transform after application of a Hanning window, over a sliding window of 512 ms. To improve reliability, P(ω, t) is taken as the average of three successive windows with 90% overlap. We used ω_{0} = 10 Hz, ω_{max} = 100 Hz, ω_{a} = 10 Hz, and ω_{b} = 30 Hz. This is similar to the SNR criterion introduced in Hurtado et al. (2005) for basal ganglia oscillations in the tremor frequency range. We set the SNR threshold value to 2 in this study. Thus SNR(t) provided us with an “instantaneous” measure for the presence of the beta-band oscillations. If the gap in time between intervals with SNR(t) above the threshold was smaller than 256 ms, it was considered to be a single continuous interval. Only intervals with SNR above the threshold were considered in the subsequent analysis.

The analysis included computation of the phase synchronization index, as described in Hurtado et al. (2004, 2005). Phase synchronization (phase locking) is a very general phenomenon in oscillating systems (Pikovsky et al. 2001) and provides an efficient tool to describe neuronal dynamics (e.g., Varela et al. 2001). As in Hurtado et al. (2004), the phase of the signals was recovered with a Hilbert transform and the synchronization index is defined as *t*_{k}. For most of the cases we considered two different window lengths 1 and 1.5 s. To estimate the significance of γ, we used surrogate data, which preserve important spectral properties of the original data, as described by Hurtado et al. (2004). We used “S4-type” surrogates, which incorporate the phase slips (Hurtado et al. 2004). An individual set of surrogates was generated for every analysis window. To compute γ, we took time-intervals of greater than threshold SNR, which were >2,000 ms. For the analysis of synchronization dynamics, we further considered only those episodes, during which γ exceeded the 95% confidence level, generated from surrogate data. Thus we obtained episodes of significant beta-band oscillatory activity, during which statistically significant synchronization occurs. The properties of this synchronization were studied via the analysis of first-return maps, a traditional tool in dynamical systems theory (e.g., Strogatz 1994).

### First-return maps for the phases and desynchronization events

First, we set up a check point for the phase of LFP (we used φ_{LFP} = 0 as a check point, but the specific value is not crucial for the analysis). Whenever the phase of the LFP signal crossed this check point from negative to positive values, we recorded the values of the spiking signal phase, generating a set of consecutive phase values {φ_{spikes,i}}, *i* = 1, …, *N*, where *N* is the number of such level crossings in an episode under analysis. Then we plot φ_{spikes,i+1} versus φ_{spikes,i} for *i* = 1, …, *N* − 1. The properties of these plots will characterize the dynamics of intermittent synchronization. Let us note that although we represent (φ_{spikes,i},φ_{spikes,i+1}) space in figures as a rectangle, this is in fact a torus because the phase is defined modulo 2π. A fully synchronous regime would result in a very simple first-return map—a single point on the diagonal φ_{spikes,i+1} = φ_{spikes,i}. Completely uncorrelated phases of the signals would yield a (φ_{spikes,i,}φ_{spikes,i+1}) space homogeneously filled with the dots. A tendency for predominantly synchronous dynamics will appear as a cluster of points, with the center at the diagonal φ_{spikes,i+1} = φ_{spikes,i}.

Since we include only those periods of activity, during which the synchrony is present, we always had this kind of a cluster. We determine the center of the cluster (as a mode of the 10-bin histogram of the phases φ_{spikes,i}) and then shift all values of the phases, to position the center of the cluster at a point with the coordinates (π/2,π/2)—the center of the first quadrant. This allows for uniform analysis of all first-return maps. We then consider how the system leaves the synchronization cluster and its vicinity and how it returns back to synchronization by quantifying transitions between different quadrants of the (φ_{spikes,i},φ_{spikes,i+1}) space (see Fig. 1 for the schematics of this space). The standard convention is to number the quadrants in a counterclockwise direction, but the dynamics in the (φ_{spikes,i},φ_{spikes,i+1}) mostly follows clockwise pattern. Thus we refer to four regions of the phase space, which are numbered in a clockwise manner. While the system is in the first region, we consider it to be in a synchronized state. It spends most of the time there. Dynamics outside of the first region will be called *desynchronization event*. We look at the following “rates”: *r*_{1} is the ratio of the number of trajectories escaping the first region for the second region (*N*_{1} → *N*_{2}) to the number of all points in the first region (*N*_{1}); *r*_{2} is the ratio of the number of trajectories escaping the second region for the fourth region (*N*_{2} → *N*_{4}) to the number of all points in the second region (*N*_{2}); *r*_{3} is the ratio of the number of trajectories escaping the third region for the fourth region (*N*_{3} → *N*_{4}) to the number of all points in the third region (*N*_{3}); *r*_{4} is the ratio of the number of trajectories escaping the fourth region for the first region (*N*_{4} → *N*_{1}) to the number of all points in the fourth region (*N*_{4}).

Because the phase space (φ_{spikes,i},φ_{spikes,i+1}) essentially represents current phase versus future phase, the transitions in this space cannot be completely arbitrary. For example, we cannot get from the second region to itself because the points in the second region are the points with low future phase and high current phase. The future phase is the current phase on the next cycle (this is the nature of the first-return maps), it is low and thus on the next cycle the point will be in the third or fourth regions. Similarly we can get from the fourth region (negative φ_{spikes,i}, positive φ_{spikes,i+1}) only to first and second regions (both with positive φ_{spikes,i}, which is equal to the previous value of φ_{spikes,i+1} just by construction of the first-return maps). However we consider all theoretically possible transitions. Therefore the considered rates provide a fairly complete characterization of how synchrony breaks down and emerges again in time.

All rates vary between 0 and 1. *r*_{1} tells how frequently synchronization is lost, the closer *r*_{1} is to 1, the more frequently the synchronized dynamics is interrupted, although it does not tell the duration of the desynchronization event. *r*_{2} characterizes the number of times the system will go to the fourth region, instead of going to the third region—an alternative synchronized state with a different phase, and *r*_{3} characterizes how long or short that third-region state is (similar to *r*_{1}). *r*_{4} characterizes the number of times the system will go back to the first region returning to the synchronized state: low values of *r*_{4} indicate that system frequently returns to the second region and the desynchronization events may be long. The predominance of short desynchronizing events requires large values of *r*_{3} and *r*_{4} or large values of *r*_{2} and *r*_{4}. In principle, *r*_{3} may take small values, which (together with small values of *r*_{1}) would indicate the presence of two synchronized states with different phases. This would correspond to a dynamics with two phase-locked states with different phase shifts, replacing one another during an episode. However, as we will show in results, this case was not observed in our data.

Our division of the phase space into four regions creates a coarse-grained picture of the phase space, but it allows for an efficient inspection of the fine temporal details. Time-averaged measures of synchrony (even if computed over relatively short time windows) characterize whether the synchronization is strong or weak overall. Utilization of these rates lets us explore the dynamics of synchronization in time and, in particular, permits evaluation of whether the weakness of synchrony is due to a few relatively long desynchronization events or to large number of short desynchronization events.

This choice of division of the phase space into four regions is not the only possible choice. Coarse-grained partition of the phase space into four areas essentially establishes ±π/2 brackets for what is considered to be similar phases. This is a compromise value, which is substantially smaller than the full cycle. Beyond these brackets, the phase difference will be closer to the full cycle. Selection of a much smaller interval would result in an overestimate of desynchrony in noisy biological data. Using much larger intervals would render very large phase deviations as nonessential. The fact that the durations of desynchronization events recovered from the data and from the rates are similar (see results) provides further argument in favor of division of the phase space into four parts. It is also important to note that while these ±π/2 brackets will affect the values of the rates, we know that the data subjected to the analysis is synchronous regardless of the value of these brackets. Selection of the synchronized episodes is done before the first-return plot analysis, so that the definition of ratios, partition of the phase space into four regions etc., had no effect on it.

It should be emphasized that the methods described allow for an analysis of the temporal development of synchrony, which is not perfect: that is the phase difference between “synchronous” signals may fluctuate. However, these fluctuations are moderate so that, overall, the selected episodes have some temporal coordination. This analysis of the fine temporal structure of the synchrony would make no sense, if there were no synchrony between the signals at all.

## RESULTS

### Data used for analysis of intermittently synchronous dynamics

The patient's recordings were processed as described in methods and yielded the following number of episodes included in the subsequent analysis: 32 episodes of average duration 4.4 ± 5.9 s if the window length for the computation of synchronization index γ was 1 s and 25 episodes of average duration 5 ± 6.7 s if this window length was 1.5 s. Although these two sets of selected data are similar, they are not the same and even larger differences may occur for larger differences in the γ-computation window length. The fact that the duration of this window affects the inclusion of the data to a certain degree is not surprising. The length of this window defines the time scales over which γ detects synchrony on average (truly instantaneous synchrony is impossible to define), and there is no rigorous way to determine the best window length. However, qualitatively, the results described in the following text do not depend much on the length of the window, and we show this by reporting the results for both window lengths. Given the relatively small number of good episodes per patient available, we pool all episodes together for the analysis. Note that overall we have a relatively short duration of good beta-band activity, which may be attributed to the relatively short intraoperative recordings and specific selection criteria. However, the purpose of this study is to investigate the fine temporal properties of synchronization of this activity not the temporal variations of the activity itself.

### Examples of the observed neuronal dynamics

To illustrate synchronous dynamics, first, we would like to present an example with a very short piece of the real data. Figure 2 illustrates both unfiltered and filtered spikes (*A*) and LFPs (*C*). The middle plot (Fig. 2*B*) represents the sine of the phase of both oscillations. Note that the amplitude of both waveforms is equal to one, so that the temporal relationship between these signals is representative of the phase synchronization. Black stars indicate the phase of one signal, as the other signal goes through the check point (see *First-return maps for the phases and desynchronization events*). Figure 3 is the first return map obtained from the data in Fig. 2. The signals in the middle plot of the Fig. 2 are not in perfect synchrony but are clearly coordinated in the time domain. As a result, all the points in Fig. 3 are positions in the first region (see *First-return maps for the phases and desynchronization events*).

Now let us consider a more complicated case, which is the subject of the study. An example of a first-return map of the phases for a long episode of oscillatory activity is presented in Fig. 4. There is a clear tendency of points to form a large cluster with the center near the diagonal. This reflects an overall tendency for synchrony, as it should, because we select episodes where some synchrony is present (see *Data analysis*). The cluster has a relatively large size, which means the synchrony is not perfect. In other words, the phase difference between the phases of spikes and LFP does not vary much with time, yet it is not perfectly constant and fluctuates in some interval (due to noise, dynamical effect, short-term plasticity, synaptic inputs from other parts of the nervous system or other factors). The same should be true for the underlying cellular and synaptic activity, which gives rise to extracellular spikes and LFP. To show the dynamics in the phase space, we draw a line from each point in all four regions of the phase space to the next successor point. Thus we see that in this example, for most of the time, a point in the second region transitions to the fourth region and from the fourth region to the first region (that is, back to synchronized state). Such a situation indicates the predominance of short desynchronizing episodes, as the system does not wander much between regions and returns to the synchronous state quickly. We report the quantitative results in the next subsection.

### Transition rates between different parts of the phase space

To characterize the temporal dynamics of the oscillatory synchronization we quantify the transitions between different parts of the first-return map. The rates *r*_{1}, *r*_{2}, *r*_{3}, *r*_{4} (see *First-return maps for the phases and desynchronization events*) describe the transition rates between different parts of the first-return plot and thus characterize the chances of transition of a system from one state to another. The values of the transition rates, averaged over all acceptable episodes of activity are presented in Fig. 5. Transition rate *r*_{1} is statistically different from other transition rates *r*_{2,3,4}, as revealed by Mann-Whitney *U* test, yielding (for 1 s long time window used to compute the synchronization index γ) significance levels *P* = 2.52*10^{−9} (*r*_{1} vs. *r*_{2}), 1.03*10^{−8} (*r*_{1} vs. *r*_{3}), and 2.88*10^{−8} (*r*_{1} vs. *r*_{4}).

The values of the transition rates *r* in Fig. 5 were computed from episodes of varying length for each episode and then averaged over all episodes. Thus the short episodes have the same impact on the mean values of *r* as long ones. To obtain transition rate average values with the long episodes having proportionally larger impact, we computed weighted average values with the weight of an episode equal to the number of points in the first-return plot for the episode (that is, approximately proportional to the duration of the episode). These results are presented in Fig. 6. The weighted transition rate *r*_{1} is also statistically different from other weighted transition rates *r*_{2,3,4}, as revealed by Mann-Whitney *U* test, yielding (for 1 s long time window used to compute the synchronization index γ) significance levels *P* = 8.08*10^{−6} (*r*_{1} vs. *r*_{2}), 1.38*10^{−5} (*r*_{1} vs. *r*_{3}), and 1.59*10^{−4} (*r*_{1} vs. *r*_{4}).

Comparison of Figs. 5 and 6 indicates that the transition rates *r*_{1} do not depend strongly on the method of computation. In principle, different approaches to averaging are not necessarily expected to produce statistically identical results. However, overall, the transition rate *r*_{1} is of the order of 0.3–0.4 and other rates are of the order of 0.6–0.7, and these values do not depend much either on the inclusion criteria (the length of the time-window used to compute the synchronization index γ) or on the method of averaging. Thus we see that on average, synchronous dynamics in STN in our patient group is interrupted every three periods of oscillations. On the other hand, the higher values of other transition rates indicate that the return to the synchronized state is relatively quick. We show below that values of these rates in the range of 0.6–0.7 translate into the most frequent desynchronization event being of very short duration. These may be characteristic features of beta-band oscillations in parkinsonian STN. The quantitative estimation of the return to the synchronized state is considered in the next subsection.

### Duration of desynchronization events

Although the transition rates presented in the preceding text provide us with the chances of transition from one dynamical state to another, they do not directly represent the duration of the desynchronized event (which may consist of a sequence of many nonsynchronous states and synchronous states with phase lags different from the major synchronous state). Given the first-return map approach we employ, the duration of the desynchronization event is the number of cycles of oscillations during which the signals are not synchronized (as judged by the time instant when one signal passes through a check point). There are two ways to estimate the duration of the desynchronized event. One is to calculate it directly from the data, computing the relative frequencies of desynchronization events of different length (which start and end at the synchronized state). The other approach is to estimate it using the transition rates: assume that the transitions between the four parts of the phase space are independent, and then each path from the second region (where desynchronization begins) to the first region (primary synchronous state) has a probability equal to the product of the transition rates for that path. For example, to estimate the probability of a desynchronization event of the shortest duration we should consider the shortest path 2-4-1, the corresponding probability is *r*_{2}·*r*_{4} (note that desynchronization will always start at the 2nd region by the virtue of our description of the dynamics). This will correspond to the desynchronization length of one cycle (with some possible rare exceptions of unusually slow or fast phase advances in one oscillatory signal with respect to the other). If we want to consider the chances of the duration of three cycles, we should consider two possible paths: 2-4-2-4-1 and 2-3-3-4-1. The probability of a desynchronization event of the duration of three cycles will be equal to the sum of products of the corresponding transition rates *r*_{3}·(1 − *r*_{4})·*r*_{2}·*r*_{4} + (1− *r*_{2})·(1 − *r*_{3})·*r*_{3}·*r*_{4}.

These computations can be easily done for the paths of short duration. The results of the estimation of the duration of desynchronization events and their direct calculation from the data are presented in Fig. 7. The latter was done in two different ways to allow for comparison with estimations based on unweighted and weighted ratios. First we computed the frequencies of desynchronization events of different duration for each episode and averaged the results across episodes. This is equivalent to unweighted rates. Second we computed the frequencies of desynchronization events of different duration for all episodes together. This is equivalent to weighted rates (each episode makes an impact roughly proportional to its length). The difference between the frequency of occurrence of the shortest desynchronization and the frequency of occurrence of the longer desynchronization events is statistically significant. The Mann-Whitney *U* test yields the significance levels *P* = 1.37*10^{−7} (duration 1 vs. duration 2), 1.12*10^{−7} (duration 1 vs. duration 3), 1.35*10^{−7} (duration 1 vs. duration 4), and 1.36*10^{−8} (duration 1 vs. duration 5) for unweighted rates from the data and *P* = 2.14*10^{−6} (duration 1 vs. duration 2), 7.41*10^{−6} (duration 1 vs. duration 3), 7.72*10^{−8} (duration 1 vs. duration 4), and 1.11*10^{−7} (duration 1 vs. duration 5) for weighted rates from the data.

Thus the results at the Fig. 7 confirm that the short desynchronization events dominate the dynamics. They are more frequent than other events by a factor of 2–3. For each of the methods, the probability of particular duration is an almost monotonically decreasing function of the duration. Qualitatively, the results do not depend on the kind of averaging procedure used and are similar for direct calculation of the frequencies of different durations and estimates based on the transition rates. This indicates that considering the rates of transitions between different regions of the first-return plot to be independent provides an adequate description of synchronization-desynchronization dynamics. Nevertheless the ratios-based estimate appears to underestimate the chances of the shortest desynchronization events and the chances of long desynchronization events. This may indicate that transitions between different parts of the phase space may depend weakly on the past dynamics. However, this effect is small. It appears that the qualitative nature of the results is more assuring, while less emphasis should be placed at particular numerical values, which may be slightly different for other patient samples and other values of parameters in the data analysis procedures.

## DISCUSSION

### Detection and temporal structure of the beta-band synchronous oscillations in STN

Our results reveal a complicated temporal structure of synchronization of beta-band oscillations in the subthalamic nucleus of parkinsonian patients. Synchronization is not an instantaneous phenomenon (see, e.g., Pikovsky et al. 2001). In the case of many real-world systems, it may exist on average, in a statistical sense. We may detect the intervals of synchrony in our data with relatively high (of the order of 1 s) resolution in a statistically rigorous way. However the dynamics may pass in and out of synchrony very briefly following a certain temporal structure. This is what we observed in our data.

The inclusion criteria allowed us to concentrate on the episodes containing synchronous activity, however the synchronization index γ does not provide information about the fine temporal structure of synchrony because the analysis window length must be sufficiently large. This window size must be long enough to provide powerful statistics and yet short enough to get reasonable time resolution. A window size of 1 s corresponds to ∼20 cycles of beta-band oscillations. Exploration of synchrony patterns on a finer time scale necessitated the construction of a somewhat unconventional time-series analysis approach: we constructed first-return plots of the phase difference between signals and studied the dynamics of the phase difference between extracellularly recorded units and LFP in this phase space.

It turned out that synchronized dynamics is interrupted by desynchronization events. These events are apparently irregular, although not completely random—there is a predominance of short desynchronization events. The signals go out of phase for just one cycle of oscillations more often than for two or a larger number of cycles. An alternative scenario (desynchronization events are longer but less frequent) would produce the same degree of average synchronization. However, this study shows that this alternative is not realized in the parkinsonian BG. Although the desynchronization event here is considered relative to a particular data analysis method, the presence of short, frequent episodes of diminished or absent synchrony should persist regardless of the details of the time-series analytic procedures. The predominance of short yet relatively frequent desynchronization events may have certain implications for the functional aspects of BG physiology in PD, which we discuss later in this section.

The quantities (transition rates), which characterize the dynamics of synchronized/desynchronized events depend on certain parameters of the time-series analysis techniques. We tested different inclusion criteria and different ways of computing the transition rates. There is no one way to select the best time-series analysis procedure. However, the resulting values of the rates and hence the fine temporal structure of intermittent synchronization do not depend much on these details. Thus our results appear to be robust and to provide a general characterization of BG dynamics in PD.

The signals under consideration are extracellular single units and LFP. Extracellular units represent somatic/axonal activity, in other words—output signal. LFP are primarily generated by synaptic potentials and thus reflect incoming and local processing activity (Buzsaki et al. 2003; Mitzdorf 1985). Similarly to the cortex, LFP recorded in different parts of BG are of synaptic origin and are relatively locally generated (Brown and Williams 2005; Goldberg et al. 2004; Magill et al. 2004). LFP are not as local as the extracellular spikes signal (Goldberg et al. 2004), but they are generated in a small vicinity of the recording electrode. The relationship between spikes and LFP in STN reveals the relationship between synaptic/dendritic and somatic/axonal activity in STN. The existence of local connections within STN is very unlikely (Wilson et al. 2004). Thus the correlations between LFP and spikes in STN may reflect the relationship between the dynamics of STN inputs and outputs.

### Mechanisms of the intermittent synchronous oscillations in BG

Synchronized beta-band oscillatory activity in STN and related structures is widespread, in general suppresses movement, and is suppressed by dopamine (see Introduction). Beyond subthalamus and pallidum functionally significant oscillations have been observed in cortico-striatal networks (Berke et al. 2004; Courtemanche et al. 2003), and the lack of dopamine promotes their strength and synchrony (Costa et al. 2006). The localization of discrete BG oscillators is questionable (Montgomery 2007), and the oscillations may arise in very distributed networks. However, the dynamics of oscillations in different parts of BG in response to dopamine and movement seems to be very consistent, so the results of the study of synchronous oscillation in one node of BG-thalamocortical networks (such as STN) may be quite typical for BG networks as a whole.

Note that the casual link between beta-band oscillations and movement deficit does not fit with the results of some studies of primate and rodent models of PD (Leblois et al. 2007; Mallet et al. 2008b), which started to observe synchronized oscillatory activity *after* detection of motor impairment. However, motor deficits in the aforementioned experiments could be of dystonic, not Parkinsonian nature (Brown 2007; Mallet et al. 2008b) and the detection of variable weakly synchronized oscillations is difficult (Leblois et al. 2007). The intermittent nature of synchronous oscillations may provide a further explanation for detection difficulty.

What may induce intermittency of synchronization of neuronal oscillations? Several factors can potentially promote the transient nature of synchronous oscillatory dynamics. A difference in frequencies of oscillators may induce intermittent synchronization (e.g., Ermentrout and Rinzel 1984). Short-term plasticity, which was observed in different parts of BG (Di Filippo et al. 2009; Fitzpatrick et al. 2001; Hanson and Jaeger 2002; Mahon et al. 2004), and fluctuating sensory input may also be desynchronizing influences. Variability of synchronization can also arise intrinsically in oscillatory systems due to moderate strength of coupling. Further studies are needed to explore the mechanisms of the observed intermittency.

There are also some experimental aspects, which may potentially contribute to the variability of time-sensitive measures of synchrony, such as patients' movement, main line (60 Hz) artifacts, cardiac and respiratory artifacts, or other infidelities of recording. However, this is unlikely to be the case here. The episodes included in the analysis were recorded when the patients were at rest. 60 Hz and very low frequency (cardiac/respiratory etc.) artifacts are spectrally far away from the beta-band and are reliably filtered out. The use of phase synchronization (which is not sensitive to the amplitude of the signals, and which, in our case, is not very sensitive to the occurrence of spikes within beta-band burst) further contributes to the robustness to the results. The use of dexmedetomidine during surgical procedure may also have some unknown effect on the observed phenomenon.

Finally, it is interesting to note, that although another type of oscillations in BG in PD—tremor—may be fundamentally different from the beta oscillations in BG (Rivlin-Etzion et al. 2008; Weinberger et al. 2006), the synchronization properties of tremor oscillations are also intermittent (Hurtado et al. 2005, 2006).

### Potential functional aspects of intermittent synchrony

The basic feedback circuits (Bevan et al. 2002; Plenz and Kitai 1999; Terman et al. 2002) and rich membrane properties of BG neurons (Bevan et al. 2006; Surmeier et al. 2005) can facilitate oscillations. BG-thalamocortical loops may be able to generate and facilitate oscillations too. Recent studies focus on the potential cellular and network mechanisms of oscillations in STN (Mallet et al. 2008a). However, very robust, perfectly synchronous oscillations may be hard to modulate (in comparison with weakly synchronous dynamics). If this is the case, they will be less efficient in transmitting information and therefore are more likely to be pathological. As we noted in the preceding text, the dynamics of STN in the normal state shows much less oscillatory synchrony. Yet even in PD the synchrony is intermittent. This intermittent, nonpermanent character may be what permits the diseased BG to perform some informational functions.

There are two primary ways to produce synchronized oscillations: coupling between oscillators and common input. Both are possible in the parkinsonian BG. The low-dopamine state is characterized by the loss of specificity of neuronal responses (Boraud et al. 2000; Filion et al. 1988), which may result in an increase in the common input into different BG oscillators. The low-dopamine state also, in general, increases the strength of synaptic transmission in BG because dopamine tends to suppress many types of BG synapses. For example, dopamine inhibits GABA release in STN (Bauferton and Bevan 2008; Cragg et al. 2004; Shen and Johnson 2005) and GPe (Cooper and Stanford 2001; Ogura and Kita 2000) and suppresses STN to GPe transmission (Hernandez et al. 2006). Thus the coupling strength becomes effectively stronger in the parkinsonian state. Yet this strong coupling and common input is not too strong as it is unable to induce complete synchrony. The intermittent synchrony that we observe in the parkinsonian case may be a result of a propensity of BG circuits to be engaged in the brief synchronized episodes of activity needed for movement control. This transient character of dynamics may be typical for neural systems (Rabinovich et al. 2008). The low-DA state with stronger coupling and stronger common input may result in a partial suppression of this very transient (and hard-to-detect) character of neuronal dynamics, favoring only short desynchronization events, which interrupt mostly synchronous episodes.

Finally, suppression of beta-band synchronization may be a potential way to treat PD motor symptoms. Thus it is important to know the dominance of desynchronization events on short time scales. It may indicate that even though the synchronization in parkinsonian BG is fragile enough to be frequently destabilized, it nevertheless can reestablish itself very quickly. This kind of knowledge may assist attempts to realize adaptive, close-loop deep brain stimulation (DBS), such as the recent suggestions for controlling synchronous oscillations by delayed nonlinear feedback (e.g., Popovich et al. 2006; Tukhlina et al. 2007).

## GRANTS

This study was supported by Indiana University Research Support Funds Grant and National Institutes of Neurological Disorders and Stroke Grant 1R01NS-067200 (NSF/NIH CRCNS program).

## ACKNOWLEDGMENTS

We thank Dr. Steven Schiff for help with the data analysis and Drs. Thomas Witt and Richard Friedman for support of the data collection.

- Copyright © 2010 the American Physiological Society