## Abstract

Most neurons in the external and internal segments of the globus pallidus and the substantia nigra pars reticulata (GPe, GPi, and SNr) are characterized by a high-frequency discharge (HFD) rate (50–80 Hz) that, in most GPe neurons, is also interrupted by pauses. Almost all (∼90%) of the synaptic inputs to these HFD neurons are GABAergic and inhibitory. Nevertheless, their responses to behavioral events are usually dominated by increases in discharge rate. Additionally, there are no reports of prolonged bursts in the spontaneous activity of these cells that could reflect their disinhibition by GPe pauses. We recorded the spontaneous activity of 385 GPe, GPi, and SNr HFD neurons during a quiet-wakeful state from two monkeys. We developed three complementary methods to quantify the balance of increases and decreases in the spontaneous discharge of HFD neurons and validated them by simulations. Unlike the behavioral evoked responses, the spontaneous activity of pallidal and SNr neurons is not dominated by increases. Moreover, the activity of basal ganglia neurons does not include bursts that could reflect disinhibition by the spontaneous pauses of GPe neurons. These findings suggest that the discharge increase/decrease balance during a quiet-wakeful state better reflects the inhibitory input of the HFD basal ganglia neurons than during responses to behavioral events; however, the GPe pauses are not echoed by comparable bursts either in the GPe or in the output nuclei. Changes in the excitatory drive of these structures (e.g., during behavioral activity) thus may lead to a remarkable change in this balance.

## INTRODUCTION

Many neurons of the basal ganglia are characterized by their spontaneous high-frequency discharge (HFD) rate. These include the majority (85%; DeLong 1971) of GPe (external segment of globus pallidus) neurons, which play a central role in the internal stages of the basal ganglia circuitry (Arkadir et al. 2004; Bolam et al. 2000; Kita 2007) and almost all GPi (internal segment of globus pallidus) and SNr (substantia nigra pars reticulata) neurons, which compose the output stage of the basal ganglia network (Gerfen 2004; Haber and Gdowski 2004; Schultz 1986). The high-frequency discharge of many of the GPe neurons is also interrupted by long intervals of total silence or pauses (DeLong 1971; Elias et al. 2007).

The long dendrites of neurons in the pallidal complex (GPe and GPi) are densely covered with synaptic boutons (Francois et al. 1984). The majority of these (>80%) represent striatopallidal GABAergic terminals and only 5–10% are subthalamic nucleus (STN) glutamatergic terminals (Difiglia et al. 1982; Shink and Smith 1995). In contrast to this anatomical dominance of inhibitory input, physiological studies consistently report that most (60–80%) of the pallidal neurons increase their firing rate in response to behavioral events (Brotchie et al. 1991; Georgopoulos et al. 1983; Jaeger et al. 1995; Mink and Thach 1991; Mitchell et al. 1987; Turner and Anderson 1997). Neurons of the SNr belong to the same morphological type as pallidal neurons (Fox et al. 1974; Yelnik et al. 1987), receive similar anatomical inputs (Francois et al. 1987; Rinvik and Grofova 1970), and as in the pallidum their behavioral-related activity is not dominated by a decrease in firing rate (Handel and Glimcher 1999; Nevet et al. 2007; Sato and Hikosaka 2002; Schultz 1986). Additionally, the soma of neurons in the GPi and the SNr (the output structures) as well as the GPe itself is innervated by GPe GABAergic synapses (Kincaid et al. 1991; Sadek et al. 2007; Smith and Bolam 1989). However, neurons in the output structures do not exhibit dominance of discharge decreases in response to behavioral events that counter the activity in the GPe (Anderson and Horak 1985; Mitchell et al. 1987; Turner and Anderson 1997).

Recent pharmacological and physiological studies appear to show that the spontaneous activity of pallidal cells is under the influence of inhibitory and excitatory inputs (Kita et al. 2004; Tachibana et al. 2008). However, the ratio between these opposing forces (“the discharge increase/decrease balance”) has not been addressed in pallidal and SNr cells during spontaneous activity, either before or after depletion of dopamine by the MPTP neurotoxin. The present study therefore aims *1*) to characterize the discharge increase/decrease balance in spontaneous activity in the GPe, GPi, and the SNr; *2*) to examine whether basal ganglia HFD neurons exhibit discharge increases that could reflect disinhibition by pauses of GPe neurons; and *3*) to compare the discharge increase/decrease balance in the normal and the MPTP states. We therefore developed three complementary methods that were validated by extensive simulations of spike trains. We used these methods to analyze the neuronal activity of HFD neurons in the basal ganglia of two monkeys in the normal and MPTP states.

## METHODS

### Recording procedures

The experiments were carried out on two monkeys (monkey Cu: *Cercopithecus aethiops aethiops*, female, weighing 3.8 Kg; monkey P: *Macaca fascicularis*, female, weighing 3 Kg). Monkey Cu was trained to perform a self-initiated button-pressing task to obtain a liquid reward. However, most of the recordings in monkey Cu were conducted during a “quiet-wakeful” state and only nonbehaving periods were included in our database (see following text). Monkey P was not engaged in any behavioral task and was trained only to sit quietly in the primate chair.

After training, a recording chamber was attached to the monkeys' skulls. The recording chamber was tilted laterally in the coronal plane and was positioned by a stereotaxic device to cover most of the pallidal (GPe and GPi) and the SNr territories (Contreras et al. 1981; Martin and Bowden 2000; Szabo and Cowan 1984). The exact position of the chamber was established using a magnetic resonance imaging (MRI) scan and electrophysiological mapping. Details of the surgery and data recording methods are given elsewhere (Elias et al. 2007; Heimer et al. 2006). Neurons in all structures were identified according to their stereotaxic coordinates (based on the MRI and the primate atlas data) and their real-time physiological identification (for example, typical spike shape and firing characteristics of neighboring neurons; for details see Elias et al. 2007).

The monkeys' care and surgical procedures were in accordance with the *National Institutes of Health Guide for the Care and Use of Laboratory Animals* (1996) and the Hebrew University Guidelines for the Use and Care of Laboratory Animals in Research, supervised by the Institutional Committee for Animal Care and Use.

### MPTP injections

Both monkeys were systemically treated with 1-methyl-4-phenyl-1,2,3,6-tetra hydropyridine (MPTP) after the recordings in the normal state. The treatment course included five intramuscular injections of 0.4 mg/kg of the MPTP neurotoxin (Sigma, Rehovot, Israel) over a period of 4 days. In both monkeys, severe Parkinsonism developed within 5 days from initiation of the MPTP injections, and recordings were resumed 4 days (monkey Cu) or 5 days (monkey P) after the last injection. The initial symptoms of monkey Cu included bradykinesia, flexed posture, and lower limb dystonia, which evolved into akinesia and freezing (days 3–4 of the MPTP treatment). Tremor and rigidity subsequently appeared (day 5 of the MPTP treatment). The initial symptoms of monkey P included flexed posture and limb dystonia (days 2–3 of the MPTP treatment), which evolved into akinesia and rigidity. Both monkeys were clinically assessed on a regular basis using primate clinical scales. The clinical state of severe parkinsonism remained stable during all days of recording with an average parkinsonism score of 30/36 in money Cu (modified primate clinical scale) and 19.5/29 in monkey P (the Kurlan scale; Imbert et al. 2000; Kurlan et al. 1991). In monkey Cu recordings were also conducted after treatment with dopamine replacement therapy (these data are not included herein).

### Neuronal database

Our database includes only stationary (rate stability; Tuckwell 1998) and well-isolated neurons that had been recorded for ≥120 s. The isolation quality was assessed both by the refractory period (Fee et al. 1996) and by a waveform isolation score that quantifies the validity of our on-line attribution of extracellular waveforms to single cells (Joshua et al. 2007). Only units with an isolation quality ≥0.8 and fraction of short interspike intervals (ISIs, <2 ms) out of the total ISIs <1.5% were included. In addition we included only neurons with an average firing rate exceeding 20 spikes/s (“high-frequency discharge”) to exclude GPe LFD-B and SNc neurons (DeLong 1971; DeLong et al. 1983). Finally, we excluded neuronal segments that were recorded during task performance (in monkey Cu). The longest continuous time segment (a single segment per neuron) that fulfilled all the above criteria was included in the neuronal database. Mean recording duration and mean number of spikes/cell per monkey, area, and state ranged from 5 to 24 min and 18,086 to 101,310 spikes (Table 1).

GPe neurons recorded in the normal state were defined as pausers or nonpausers according to objective criteria (Elias et al. 2007). Briefly, pauses were defined by a modified surprise algorithm triggered by core ISIs >250 ms. GPe cells were defined as pausers if ≥80% of their 1-min segments of recordings contained at least two pauses. Most (79.6%, 74/93) GPe HFD neurons were defined as pausers. We analyzed pausers and nonpausers separately, as well as the whole GPe population. Neurons of the output structures (the GPi and the SNr) were analyzed separately, but were also pooled and analyzed as a single group. The data for each monkey were analyzed separately, but neurons from the same structure and the same state were also pooled and analyzed as one group. Table 1 provides a detailed description of the neuronal database used here.

### Quantifying the discharge increase/decrease ratio

Four general considerations were taken into account in analyzing the neuronal data.

Even after many years of research, there is no agreement as to the optimal method for burst detection (see, e.g., Cocatre-Zilgien and Delcomyn 1992; Gourevitch and Eggermont 2007; Kaneoke and Vitek 1996; Legendy and Salcman 1985; Tam 2002). Therefore precise detection of increases and decreases in spontaneous activity is method dependent. This emphasizes the need for testing the data by different methods and by different parameters in each method.

Most available methods were designed for detecting bursts or increases in discharge rate. However, in the current study, increases and decreases had to be detected with similar sensitivity. Since our null hypothesis holds that the discharge increase/decrease balance of the spontaneous activity is not different from that of the evoked activity (i.e., >1), we prefer methods with a slight bias toward an increase/decrease ratio >1.

We assume that there are similarities between spontaneous activity and behavior-related activity (i.e., explicitly related to behavioral events triggered or recorded by the experimenters). Therefore in cases in which an a priori guess was required, we searched for spontaneous neuronal segments with parameters similar to those observed in response to behavioral events (e.g., changes in firing rate for ∼100 ms).

The validity of a method should be assessed by applying it to simulated neurons with different ratios between increases and decreases in discharge rate (Fig. 1). The simulated database should follow similar dynamics (e.g., duration and amplitude of increases and decreases in discharge rate) to that observed in response to the behavioral events. The use of synthetic data provides a rough idea of the range of different increase/decrease scores along different simulated increase/decrease ratios.

We used three methods to quantify the discharge increase/decrease balance of our data. The first method is based on modeling the neural activity with hidden Markov models (HMMs). For each neuron we estimated the locally optimized model parameters and defined the increase/decrease score as the ratio between the cumulative probability of the increase state and the cumulative probability of the decrease state. The second method is based on comparing the total weight of segments with a firing rate higher than the median firing rate to the total weight of segments with a firing rate lower than the median firing rate. The last method, the extended surprise method, uses principles similar to the surprise burst-detection algorithm (Legendy and Salcman 1985). We extended the basic surprise method for detection of discharge decreases and for analyzing the data along a range of surprise values. Details of these methods are given in appendices A, B, and C. We validated our analysis methods by applying them to simulated neurons of two types: Poissonian cells and Markov model neurons. The simulated neuronal database was composed of neurons with a baseline firing rate of 50–55 spikes/s and the simulated Markov model neurons also included varying ratios of discharge increases and decreases (details given in appendix D).

## RESULTS

### Neuronal and simulated database

Our neuronal database consisted of 385 HFD neurons recorded from three basal ganglia structures: the GPe, the GPi, and the SNr of two monkeys. In all, 138 neurons were recorded in the normal state (93 GPe cells, 45 GPi/SNr cells) and 247 neurons were recorded after systemic MPTP treatment and during a stable and severe parkinsonian state (174 GPe cells, 73 GPi/SNr cells). Only very well isolated neurons (isolation quality >0.8; Joshua et al. 2007) with stable firing rates and an unambiguous refractory period were included in the database. Details of recording time, firing rate, and fraction of “refractory period” ISIs are given in Table 1.

We used three methods to assess the ratio between increases and decreases of discharge rate. To validate these methods we applied them first to simulated Poissonian neurons (baseline firing rate ∼50 spikes/s, with and without an absolute refractory period) and Markov model neurons. Briefly, the simulated Markov model database was composed of HFD neurons with a baseline firing rate of 50 spikes/s embedded with segments of discharge increases and decreases. We used different ratios of discharge increases and decreases in relation to three parameters: frequency, duration, and amplitude changes. In the balanced conditions the parameters of discharge increases and decreases were identical. For example, in Fig. 1*A* we used the following parameters: mean frequency, 10/min; mean duration, 300 ms; and mean amplitude change, 30 spikes/s. Figure 1*B* presents raster plots of simulated Markov model neurons with different ratios of discharge increases and decreases. For more details of the simulation database, see appendix D.

### Basal ganglia discharge increase/decrease ratio: the hidden Markov model (HMM) method

Our initial working hypothesis held that the spiking activity of a single cell is a realization of a process with three states. Each state is characterized by a distinct firing probability. These states are a baseline state, an increase state (with a firing probability higher than the firing probability in the baseline state), and a decrease state (with a firing probability lower than the firing probability of the baseline state). It was therefore natural to apply the HMM analysis to our data. We first estimated the maximum-likelihood parameters of a three-state HMM model by applying the Baum–Welch algorithm. This algorithm starts from an initial guess of the model parameters and converges to a local maximum estimation of these parameters. We used two methods of selecting the initial guess: a fixed model for all neurons and a variable model, which was the most likely of 12 possible models (see appendix A). Based on the obtained model parameters for each neuron, we calculated the posterior probabilities of each of the states. We defined the discharge increase/decrease score as the ratio of the mean probability of being in state 2 (increase) to the mean probability of being in state 3 (decrease) (e.g., discharge increase/decrease scores <1 suggest a tendency toward discharge decreases; see the pauser cell in Fig. 2).

In the first step we applied the HMM method to our simulated neurons. In the Poissonian cells without a refractory period the mean discharge increase/decrease ratio was 1.14 ± 0.08 (Fig. 3 *A*, *left*), whereas in the case of the Poissonian cells with a refractory period (Fig. 3*A*, *right*) the values were slightly higher with a mean value of 1.45 ± 0.08 (Student's *t*-test, *P* < 0.001). In the simulated Markov-model neurons (Fig. 3*B*), the discharge increase/decrease score discriminated the different configurations of the model well. However, low scores (<1) were often attributed to balanced neurons (with similar probabilities for increases and decreases).

In our physiological data the distribution of the increase/decrease score according to the HMM method was skewed toward lower values in all nuclei we analyzed, both in the normal state and the MPTP state (Fig. 4). This may indicate dominance of decreases or a balanced condition (due to the biased results of the balanced simulated neurons). However, these results were calculated with a tolerance value of 0.01. The tolerance value controls the number of steps the Baum–Welch algorithm executes for estimating the model parameters. For example, decreasing the tolerance value should lead to models with better fit to the data and smaller fluctuations around the same parameters. However, when we decreased the tolerance value of the HMM models of our neuronal database, it resulted in an elevation of the increase/decrease score, primarily due to attributing longer segments to the increase state instead of the baseline state. For example, the mean probability of being in the discharge increase state of the pauser shown in Fig. 2*A* is 0.19 with a tolerance of 0.01, but it rises to 0.58 with a tolerance value of 0.0001. The mean probability of being in the discharge decrease state of the same cell is 0.25 with a tolerance of 0.01, and it remains the same with a tolerance of 0.0001. Figure 5 shows the effect of using different tolerance values on the decoding of similar raster plots. The results of low tolerance models do not appear to coincide with our visual assessment of increases and decreases and asymptotic results were not obtained in the tolerance range tested (models with tolerance values <0.0001 could not be estimated in a reasonable computing time and were therefore not tested).

The nonconvergent results of the three-state HMM analysis led us to examine the fit of a two-state model (baseline state and decrease state) to our normal physiological data. In this case, the results (fraction of time in the baseline state, fraction of time in the decrease state, and the baseline/decrease ratio) were much more stable along different tolerance values, especially in the GPe pauser population.

In summary, although the HMM method can assign probabilities that fit our visual assessment (Fig. 2), this method seems to be dependent on the tolerance value. Decreasing the tolerance value leads mainly to a decrease in the fraction of the background activity in the three-state model. The two-state (background, decrease) HMM model seems to better represent our neuronal data, especially the GPe pausers. We conclude that the discharge rate of basal ganglia HFD neurons does not represent a three-state phenomenon and seems to be better described as a realization of a one- or two-state (background-decrease) process. Thus the HMM analysis does not support the existence of prolonged bursts that reflect disinhibition (“mirroring”) by GPe pauses in either the GPe population or the GPi/SNr population. In the next sections we calculate the discharge increase/decrease ratio using continuous models with no assumption of discrete states of the basal ganglia discharge.

### Basal ganglia discharge increase/decrease ratio: the weighted probability density function

This method is based on comparing the total weight of discharge rates above the median firing rate to the total weight of discharge rates below the median firing rate in a weighted (by the deviation from the median) probability density function. To calculate this score we multiplied a normalized firing rate probability density function (PDF) by a weight vector. This weight vector reflects the certainty that a change in firing rate is indeed an increase (in cases of firing rates above the median) or a decrease (in cases of firing rates below the median). The increase/decrease score is the ratio between the sum of the values above the median to the sum of the values below the median in the weighted function. A demonstration of this method appears in Fig. 6 and details are given in appendix B.

We first applied the weighted PDF method to our simulated database (Fig. 7). In the simulated Poissonian neurons without a refractory period (Fig. 7*A*, *left*) the values of the discharge increase/decrease score were slightly higher than the expected value of 1 (1.11 ± 0.03, mean ± SD), thus suggesting a bias toward detection of discharge increases. This bias probably reflects the difference between the Poissonian distribution and the normal distribution and is affected by the smoothing procedure. However, the bias toward detection of discharge increases in the Poissonian simulation was balanced by the addition of the refractory period to the simulated Poissonian neurons. The distribution of the increase/decrease score in the Poissonian cells with a refractory period was better centered around 1 for all smoothing widths tested (e.g., in the case of σ = 100, Fig. 7*A*, *right*: 1.02 ± 0.03, mean ± SD; compared with the Poissonian cells without a refractory period, Student's *t*-test, *P* < 0.001). The weighted PDF method differentiated the different configurations of the simulated Markov-model neurons (Fig. 7*B*). We found significant *R*^{2} values (0.82, 0.92, and 0.93) between the model parameters and the increase/decrease score in each of the three simulation configuration groups.

The next step was to apply the weighted PDF method to our neuronal database. First, we visually assessed the normalized firing rate PDFs of the neurons and found that they tended to have a normal distribution (Fig. 6*F*, *left*). However, a clear bimodal distribution was sometimes noted (Fig. 6*C*, *left*). In line with the two-state HMM model results, the bimodal distribution was found especially in the GPe pauser population of both monkeys (25 of 42 pausers, 60% in monkey Cu; 16 of 32 pausers, 50% in monkey P). We then calculated the increase/decrease score of our neuronal data (Fig. 8). In the normal state, there was a clear difference between GPe pausers and other HFD cells. The distribution of the increase/decrease score of GPe pausers was skewed toward low values (i.e., toward decreases): mean increase/decrease score = 0.73, 0.93, and 0.82 for monkey Cu, monkey P, and both. The mean score of the pooled GPe population was significantly lower than the mean score of the Poissonian cells with a refractory period as well as the mean score of the balanced simulated set (Student's *t*-test, *P* < 0.001 in each case). Excluding GPe pausers, the distribution of the increases/decrease score of the remainder of our HFD cells (in GPe and output nuclei) was concentrated around 1, although with a relatively wide range (e.g., the mean ± SD score of the pooled neurons in the output structures of both monkeys was 1.05 ± 0.2, range: 0.63–1.53). The mean score of the pooled output neurons was smaller than the mean score of the Poissonian cells with a refractory period as well as the balanced simulated set, but significant only in the latter case (Student's *t*-test, *P* < 0.05). In the MPTP state, a change in the score toward higher values (e.g., more increases) was seen in GPe neurons of monkey P (mean score in normal state: 0.94 ± 0.21, vs. mean score in MPTP state: 1.1 ± 0.17, Student's *t*-test, *P* < 0.001). A similar although nonsignificant trend was observed in GPe cells of monkey Cu (mean score in normal state: 0.77 ± 0.2, vs. mean score in the MPTP state: 0.82 ± 0.23, Student's *t*-test, *P* > 0.05).

We also repeated the calculation of the discharge increase/decrease score with a nontruncated weight vector (see appendix B) as well as with different smoothing σ values (20, 50, and 200 ms). In all of these cases the same features of the population increase/decrease scores were preserved (data not shown). We also applied the same analysis by normalizing the firing rate PDF to the mean rather than to the median. In this case, the increase/decrease score distribution was better centered around 1, with slightly left-skewed distributions in GPe pausers (e.g., mean increase/decrease score of the pooled GPe pausers was 0.97 ± 0.05, range: 0.77–1). As with the increase/decrease ratio calculated by the median weighted PDF, increases of discharge rate in the spontaneous activity of HFD neurons did not dominate (data not shown).

We conclude that the weighted PDF method is insensitive to the extensive range of tested parameters and therefore yields a reliable estimate of the increase/decrease ratio. Using this method, we can divide the BG HFD neurons into GPe pausers and all other nonpauser cells. In both groups, the increase/decrease ratio was not significantly >1, indicating that unlike the case of behavior-related activity, the spontaneous activity of these neurons is not dominated by increases in discharge rate. However, the weighted PDF method assumes a linear weight to deviations of discharge rate from the expected (median or mean) discharge rate. We therefore applied a different method to assess the increase/decrease ratio that is based on deviations from the expected values of a Poissonian process, as discussed in the following text.

### Basal ganglia discharge increase/decrease ratio: the extended surprise method

We developed (see details in appendix C) an extended version of the Poisson surprise method (Legendy and Salcman 1985) for detection of positive and negative modulations in discharge rate. The method detects segments with firing rate changes (increases or decreases) where the probability of obtaining them under a Poissonian assumption is low. Applying this algorithm to segments of our neuronal recordings yielded a fine-grained detection (compared with our visual assessment) of increases and decreases (Fig. 9, *A* and *C*). The discharge decreases of the GPe pauser in Fig. 9*A* are characterized by high surprise values (>30). These values reflect the low probability, according to the Poissonian assumption, of pause segments containing a very small number of spikes for long durations. However, similar (mirror) segments with discharge increases (with similar surprise values) were not found in the same raster plot of the pauser cells (Fig. 9*A*) or in the raster plot of the nonpauser cells (Fig. 9*C*).

To enable appropriate comparison between discharge increases and decreases we used the same surprise threshold for both of these cases. Instead of selecting a single surprise threshold as is commonly done (see, e.g., Legendy and Salcman 1985; Levy et al. 2002; Wichmann and Soares 2006), we conducted the comparison along the entire range of surprise values of a given neuron. In each neuron we calculated the ratio between the frequency of increases and the frequency of decreases for each of the surprise thresholds. This analysis yielded a discharge increase/decrease ratio curve that was calculated for different surprise thresholds and therefore characterized the increase/decrease balance of a single neuron (Fig. 9, *B* and *D*). For example, the two neurons in Fig. 9 have an increase/decrease ratio >1 (more increases than decreases) for low surprise values; however, for most of the tested range of surprise values these two neurons display increase/decrease values of <1 (more decreases than increases). The discharge increase/decrease ratio curves for the single neurons were averaged with the rest of the neurons in the same group (e.g., cells that were recorded from the same nucleus in the same state or simulated cells with the same parameters).

To validate this method, we applied it to our simulated data (Fig. 10). In the Poissonian cells (Fig. 10*A*) the frequency of discharge increases was higher than the frequency of decreases for the different surprise thresholds. In the cells without a refractory period (Fig. 10*A*, *left*), the discharge increase/decrease ratio exhibited high values as the surprise threshold increased. This effect is probably related to the nonsymmetric nature of the Poissonian distribution: higher surprise thresholds reflect lower probabilities and in the Poissonian/binomial distribution the lowest probabilities are attributed to values above the mean (increases) rather than to values below the mean (decreases). In the simulated cells with a refractory period (Fig. 10*A*, *right*) the discharge increase/decrease ratio shifted to values <1 for surprise thresholds >3; however, the maximal surprise value was <10. We also applied this method to the simulated Markov-model neurons (Fig. 10*B*). The increase/decrease ratio curve differentiated the different configurations of the model, but usually only above a certain surprise threshold (∼8; better differentiation was noted in the frequency group; Fig. 10*B*, *inset* of *top plot*). Thus the discharge increase/decrease ratio of the Markov model simulated neurons calculated with the extended surprise method reflects their increase/decrease balance. Nevertheless, the method tends to be biased toward values >1 (i.e., biased toward detection of increases and therefore toward our null hypothesis that spontaneous increases in discharge rate are more common than decreases). Real data results, as presented in the examples in Fig. 9, with a dominance of low increase/decrease ratio values (e.g., more decreases than increases) can therefore be considered as a strong indication for dominant inhibitory control of the basal ganglia spontaneous discharge rate.

We applied the extended surprise analysis to our neuronal database (Fig. 11 *A*). In the normal state, the frequency of decreases in the discharge rate was higher than the frequency of increases for almost all the surprise thresholds. This feature was observed in all nuclei we analyzed (GPe, GPi, and SNr), but was accentuated in GPe cells. In addition, the range of the surprise values in the basal ganglia neurons (especially in GPe pausers) was much larger than the range of these values in the simulated cells. For example, a surprise value of ≥20 was found in 8.46% of the firing rate changes in GPe cells of both monkeys, but in only 0.05% of the firing rate changes in the simulated Markov neurons in the balanced condition, and not even once in the Poissonian cells with or without a refractory period. We compared the distribution of high surprise values (in which the increase/decrease ratio favors decreases in real neurons) of real neurons (pooled GPe neurons or pooled output neurons) to simulated neurons (Poissonian cells with a refractory period or balanced Markov neurons). As a statistic we used the mean of surprise values >20 and assessed statistical significance by a permutation test (*n* = 1,000). In all four comparison (real neurons vs. simulated neurons), real neurons had higher surprise values (*P* < 0.01).

In the MPTP condition, increases were more frequent than those in the normal condition. This change, which was weak in GPe cells, was clearly seen both in the GPi and the SNr (i.e., in the output nuclei of the basal ganglia). To better grasp the discharge increase/decrease balance we plotted (Fig. 11*B*) the results separately with a rigid surprise threshold (*S* = 10). This threshold is commonly used for burst detection in electrophysiological studies (Legendy and Salcman 1985; Wichmann et al. 1999). At this threshold the difference in the increase/decrease ratio between the normal and MPTP states was clear for both the GPe and the output neurons of the pooled population of both monkeys. Similar results were obtained with a surprise threshold of 15 (data not shown).

## DISCUSSION

In the current study we analyzed the spontaneous activity of 385 well-isolated HFD neurons from several basal ganglia structures: the GPe, the GPi, and the SNr during normal and MPTP states. We used three complementary methods to quantify the ratio between increases and decreases of the spontaneous discharge rate in these structures. Our main findings are that *1*) increases do not constitute the main firing rate change in the spontaneous activity of pallidal and SNr neurons, by contrast to their response to behavioral events; *2*) the firing pattern of GPe pauser cells and other BG HFD neurons does not include discharge increases that could reflect disinhibition by GPe spontaneous pauses (“antipauses”); and *3*) according to the surprise method, the activity in the output nuclei in the MPTP state is characterized by a rise in the discharge increase/decrease ratio such that increases are more frequent than in the normal state.

### Can we quantify the relations between increases and decreases of spontaneous activity?

We used three different methods to quantify the relations between discharge increases and decreases of basal ganglia “spontaneous activity.” We validated these methods by applying them to simulated cells with varying ratios of discharge increases and decreases. All the methods were able to differentiate clearly between different model configurations. However, we consider the results of two of these methods (the weighted probability density function method and the extended surprise method) as more reliable for quantifying the discharge increase/decrease ratio. In each of these two methods the principle of assigning equal weight to increases and decreases is interpreted differently: in the weighted PDF method, segments with the same deviation from the median firing rate (or the mean firing rate) are considered equivalent, whereas in the extended surprise method segments with similar probabilities under a Poissonian distribution are considered equivalent. However, neither method yielded dominant increases versus decreases in the spontaneous discharge rate. The slight bias toward detection of increases in the simulated data, especially in the surprise method, further corroborates this finding. The HMM analysis supports this view by indicating that a two-state (background and decrease states) model represents the activity of GPe pausers better than a three-state model, which also includes an increase state.

### Is the discharge increase/decrease balance during spontaneous activity similar to that observed during behavioral activity?

Anatomical and electrophysiological studies of BG HFD cells provide different estimations of their afferent excitatory/inhibitory balance. Anatomically, the majority of the afferent synapses to the pallidum and the SNr are inhibitory GABAergic synapses (Difiglia et al. 1982; Francois et al. 1987; Shink and Smith 1995). However, most of the responses of primate pallidal cells to behavioral events are increases in firing rate (Brotchie et al. 1991; Mink and Thach 1991; Turner and Anderson 1997). With respect to SNr neurons, classical studies of their saccade-related activity ascertain that it is almost exclusively composed of decreases in discharge rate (Hikosaka and Wurtz 1983a,b). However, more recent studies of SNr physiology report varying and more balanced ratios of increases and decreases in response to saccades (Handel and Glimcher 1999; Sato and Hikosaka 2002) as well as to arm and mouth movements (Magarinos-Ascone et al. 1992, 1994; Schultz 1986) or other task-related activities (Nevet et al. 2007).

Several suggestions have been put forward to explain this discrepancy, but most have been ruled out. Although γ-aminobutyric acid (GABA) may have excitatory effects (Ben Ari 2002; Choi et al. 2008; Wagner et al. 1997), the striatopallidal GABAergic pathway was shown to be inhibitory (Chan et al. 2004; Park et al. 1982; Rav-Acha et al. 2005). Another possibility is that subthalamic excitatory afferents are located at “strategic locations” such as the soma or dendritic hot spots (Hanson et al. 2004). However, subthalamic synapses are apparently evenly distributed (as striatal synapses) on pallidal cells (Gerfen and Wilson 1996; Smith et al. 1998). Moreover, GPi and SNr cells also receive inhibitory GABAergic synapses from the GPe that are strategically located on the soma (Kincaid et al. 1991; Parent and Hazrati 1995; Smith and Bolam 1989). A different possibility is that the behavioral tasks were biased toward a certain type of behavior with incongruent influence on discharge increases and decreases, as suggested by some studies (Sato and Hikosaka 2002; Turner and Anderson 2005). However, the relative lack of neural responses with decreases in discharge rate was noted in several studies with a broad range of behavioral paradigms (Brotchie et al. 1991; Handel and Glimcher 1999; Mink and Thach 1991; Mitchell et al. 1987).

We attempted to shed light on this paradox by analyzing the increase/decrease balance of the spontaneous activity of HFD cells. Although a quiet-wakeful state does not preclude some behavioral activity, it is expected to be different from a behavioral state in which there is active enrollment in a behavioral task. This is reflected for example in the GPe pausing activity that we recently described (Elias et al. 2007). The extended surprise method indicates that pallidal and SNr spontaneous activity is clearly dominated by decreases, whereas the weighted PDF method indicates a more balanced condition with a tendency toward decreases in GPe pausers. Thus according to both methods, increases are not the dominant change in the cell's firing rate. The increase/decrease balance of the spontaneous activity of BG HFD neurons is therefore different from that reported in response to behavioral events and better reflects the anatomical afferents to these cells. Our analysis does not enable us to separately address complex (e.g., biphasic) modulations in discharge that have sometimes been reported (Betarbet et al. 1997; Jaeger et al. 1995; Turner and Anderson 2005). Similarly, it is possible that the spontaneous discharge modulations are related to the oscillatory activity that has been described in the basal ganglia (Allers et al. 2002; Ruskin et al. 1999; Wichmann et al. 2002).

### Are there discharge rate increases in HFD cells that mirror the pauses of GPe neurons?

We recently analyzed the long periods of silence of GPe cells (“pauses”) (Elias et al. 2007). The somatic synaptic connections between GPe neurons (Bevan et al. 1998; Francois et al. 1984; Sato et al. 2000) raise the question of lateral disinhibition and the existence of parallel discharge increases (“antipauses”). We found no evidence for such segments in our recordings of the spontaneous activity of GPe pauser cells. Thus although periods of increases in firing rate can be found in the activity of GPe pausers (see, e.g., Fig. 9*A*), we cannot characterize the firing pattern of these neurons as a symmetric three-state model (background activity, increases, and decreases/pauses). It is still unclear whether decreases in firing rate and pauses in these cells represent a continuum of the same biological phenomenon (pauses as an extreme decrease in firing rate) or whether they are two distinct processes. The bimodal distribution of the discharge rate in many GPe pausers (50–60%, Fig. 6*C*, *left*) as well as the HMM analysis support the view that these are two distinct phenomena.

In addition, our analysis enables us to compare the activity in the GPe and the output structures (the GPi and the SNr). Classical models of the basal ganglia (Albin et al. 1989) as well as the inhibitory GPe synapses, which are located on the soma of GPi/SNr neurons (Gerfen 2004), predict reciprocal activities in the GPe and the output structures. However, our results imply that the spontaneous activity in the output structures does not mirror the activity in the GPe and particularly does not include segments of discharge increases that match the pausing activity of GPe cells. This is consistent with previous studies that report similar mean firing rates in the GPe and the GPi in the normal state as well as similar responses to striatal stimulation (Filion et al. 1991; Nambu et al. 2000; Yoshida et al. 1993).

### Can changes in the MPTP state be characterized by the discharge increase/decrease balance?

Several discharge parameters are used to compare the MPTP state to the normal state, such as the mean firing rate (Filion and Tremblay 1991; Miller and DeLong 1987), the firing pattern (bursting and oscillatory activity; Filion and Tremblay 1991; Heimer et al. 2006; Soares et al. 2004), and the degree of synchronization (Nini et al. 1995; Raz et al. 2000).

Here we used the discharge increase/decrease balance to characterize changes in the firing pattern of spontaneous neural activity in the MPTP state compared with the normal state. The results of the weighted PDF and the extended surprise are not similar in this regard, probably because they consider different aspects of the firing pattern (absolute deviation vs. Poisson surprise). We found that the discharge increase/decrease balance calculated by the surprise method in the output nuclei shifted toward higher values, such that the ratio favored decreases less. This is in line with a previous study that showed a decrease in the ratio of inhibited-to-activated GPi cells in response to passive movement after MPTP treatment (Boraud et al. 2000). This change can be explained by the hyperexcitability of the STN in MPTP state and in Parkinson's disease (Bergman et al. 1990, 1994; Vila et al. 2000).

### What is the functional significance of the spontaneous increase/decrease balance of the GP and the SNr?

Our results suggest that the transition from spontaneous activity to behavioral related activity in pallidal and SNr cells is accompanied by a major change in the discharge increase/decrease balance. This change is probably due to a powerful excitatory input from the STN, which is considered part of the fast “hyperdirect” cortex–STN–pallidal pathway. The discharge increase/decrease ratio of STN neurons during behavioral activity is estimated to be 2.4–9 (Nambu et al. 2002) and therefore can account for the change in the increase/decrease ratio of pallidal and SNr cells. In addition, we found no evidence for prolonged bursts that mirror GPe pauses. In light of the uncorrelated activity between GPe pauses (Elias et al. 2007), this result implies that a single GPe pause is not powerful enough to disinhibit its target neurons.

The functional role of the spontaneous high-frequency firing of pallidal and SNr neurons is still unclear, although several suggestions have been put forward regarding the coding advantages of this highly tonic activity (Person and Perkel 2005; Smith and Sherman 2002). Future studies of the discharge increase/decrease ratio of these structures should further contribute to our understanding of the role of the basal ganglia's spontaneous high-frequency discharge and GPe pauses in health and disease.

## APPENDIX A: THE HMM METHOD

Hidden Markov models are probabilistic models that are used in a wide range of applications (Cappe et al. 2005; Rabiner 1989), including analysis of neuronal activity (see, e.g., Branston and el Deredy 2001; Britvina and Eggermont 2007; Seidemann et al. 1996). Usually these models assume a finite set of discrete states and a set of possible outputs (observations). Each state is characterized by a distinct set of probabilities to obtain the corresponding outputs. The state is not observed directly and therefore it is “hidden.” An important property of these models (the Markov property) is that the transition probability from the state at the current time point to another state at the next time point depends solely on these two states (i.e., no memory effect).

We assume that the spiking activity of a single neuron is a realization (output) of a hidden Markov process. In any time bin (Δ*t*) the activity of a single neuron is assumed to be in one of several states (1 ≤ *i* ≤ *N*), with a probability to emit a spike *b _{i}* (and hence the probability of the second possible observation, no spike, is 1 −

*b*). The transition probability from state

_{i}*i*at time bin

*t*to state

*j*at time bin

*t*+ 1 is denoted by

*a*. The initial state probability (the probability of being in state

_{ij}*i*at the first time bin) is denoted by π

_{i}. Thus we can characterize a Markov model by the following parameters: λ = (

*A*,

*B*, π), where

*A*,

*B*, and π are the transition probability matrix (

*a*), the output probability vector (

_{ij}*b*), and the initial state probability vector (π

_{i}_{i}), respectively.

In analyzing the data we considered the activity of a single HFD neuron as a sequence of observations (outputs) with a time bin of 1 ms (Δ*t* = 1 ms). We assume that in a single time bin the neuron can be in one of three states (*N* = 3): state 1, a baseline state (with firing probability *b*_{1}); state 2, an increase state (with firing probability higher than the baseline state, *b*_{2} > *b*_{1}); state 3, a decrease state (with firing probability lower than the baseline state, *b*_{3} < *b*_{1}). We assume that in step 0 (prior to the first observation) the model is in the baseline state, and thus the initial state probabilities equal the transition probabilities from the baseline state (π_{i} = *A*_{1i}).

In the case of our physiological extracellular recording data, we only know the spike times of the neuron, but not the model parameters (λ) or the current states. There is no optimal way of estimating the model parameters (Rabiner 1989). However, we can use an iterative procedure such as the Baum–Welch algorithm for estimating the model parameters by finding a local maximum of the likelihood function. This iterative procedure, which is used for “training” the HMM, is based on an initial guess of the model parameters. After estimating the model parameters (e.g., by the Baum–Welch algorithm or simply by assuming a certain model), we are able to: *1*) calculate the likelihood of obtaining the observation sequence given the model parameters, which enables us to compare the fit of different models to a given neuron; and *2*) calculate the posterior probabilities of different states (i.e., the conditional probabilities of being in each of the states along the time units) given the model parameters (see following text).

### Quantifying the discharge increase/decrease ratio

The HMM discharge increase/decrease ratio is based on finding the optimized model parameters of each neuron and then calculating for each bin the posterior probabilities of the increase and decrease states. We defined the discharge increase/decrease score as the ratio between the cumulative increase and decrease probabilities.

The first step was to find the optimized model parameters based on a three-state model and the Baum–Welch algorithm. However, this iterative procedure is based on an initial guess of the model parameters and it converges only to a local maximum estimation. Thus different models could be obtained for different initial guesses. We therefore selected the initial guess in two ways: *1*) *Fixed guess*: for each neuron we used the same initial guess, in which the mean duration of discharge increases or decreases was 100 ms; the frequency, 10/min; and the mean amplitude change was half of the mean spike rate of the cell (e.g., for a cell that fired at 50 Hz, the initial guess of the firing rate in increase and decrease states was 75 and 25 Hz, respectively). *2*) *Variable guess*: for each neuron the initial guess model was the most likely model out of 12 models with varying ratios between discharge increases and decreases (including conditions in which discharge increases and decreases were balanced). These models were organized symmetrically so that for each model with dominant increases in relation to a certain parameter we generated a matching model with dominant decreases in relation to the same parameter. Thus we preserved the similar weight principle for detecting discharge increases and decreases. In these initial guesses we used the following parameters for discharge increases and decreases (in different combinations): frequencies of 7.5, 10, and 15/min; durations of 200, 300, and 400 ms; amplitude changes of 0.4, 0.6, and 0.8 out of the mean firing rate. To reduce the model complexity, we did not allow transitions between discharge increases and decreases or vice versa in any of the initial guess models. The results of our physiological data were virtually unaffected when the variable initial guess method was used (especially for low tolerance values), but some dependence was noted in our simulated data. However, we did not attempt to conduct an exhaustive test of different initial guesses.

After estimating the model parameters, we calculated the posterior probabilities of each of the states in each of the time bins. The discharge increase/decrease score was defined as the ratio between the mean probability of being in state 2 (increase state) to the mean probability of being in state 3 (decrease state). This score reflects the ratio between the fraction of time in increases and the fraction of time in decreases.

Our HMM analysis was designed to evaluate the increase/decrease balance of the spontaneous discharge of the neurons. It was therefore based on the assumption of a three-state model with equal weight assignment to increases and decreases. However, we also examined the effect of assuming a two-state model (a baseline state and a decrease state) on the physiological normal state data. In this case we used our “fixed initial guess” but with zero probability for transition between the baseline state and the increase state.

The HMM analysis was conducted with the statistics toolbox of Matlab 7 software (The MathWorks, Natick, MA). In the training algorithm we established two convergence criteria: *1*) the maximal number of iterations for the estimation process and *2*) the tolerance used for testing convergence of the iterative estimation process. The tolerance value controls the number of steps the Baum–Welch algorithm executes to estimate the model parameters. The algorithm terminates only when several quantities are smaller than the tolerance value. Briefly, these quantities include in Matlab software the log likelihood of the observation sequence under the current estimated model as well as the change in the transition and emission matrices between the current iteration and the previous one. We used a maximal number of iterations of 500 and different tolerance values (0.01, 0.001, and 0.0001). Unless mentioned otherwise, results are shown for the fixed initial guess and for a tolerance of 0.01 (differences in relation to other parameters are discussed in results).

## APPENDIX B: THE WEIGHTED PROBABILITY DENSITY FUNCTION METHOD

This method is based on comparing the total weight of segments with a firing rate higher than the median firing rate to the total weight of segments with a firing rate lower than the median firing rate. The contribution of each segment is weighted (by the deviation from the median) to calculate a weighted probability density function (weighted PDF). A step-by-step demonstration of this method is given for two cells in Fig. 6.

First, the firing rates of the cells were plotted as a function of time using 1-ms bins and smoothed with a Gaussian sliding window of σ = 100 ms (the Gaussian curve was truncated at 3σ = 300 ms from each side and was normalized so that the cumulative area of the Gaussian curve equaled 1). The smoothed firing rate vector was then normalized to its median value (Fig. 6, *B* and *E*) with median rate equals one. We calculated the PDF of the normalized firing rate vector in bins of 0.01 (Fig. 6, *C* and *F*, *left*). As a result of the normalization to the median, the probability of finding time segments with a smoothed firing rate >1 equaled the probability of finding time segments with a smoothed firing rate <1. Instead of setting a rigid threshold for defining increases and decreases in discharge rate (e.g., increases are segments with a smoothed normalized firing rate >1.5 and decreases are segments with a smoothed firing rate <0.5), we used a “fuzzy-logic” approach to reflect the certainty that a certain time segment was a true increase or decrease in firing rate. We multiplied the smoothed and normalized PDF with a linearly expanding probability vector ranging from 0 to 1 (the “weight vector”) in increments of 0.01 (Fig. 6, *C* and *F*, *middle*). The central bin of this weight vector, which corresponds to the central bin of the probability density function, had a value of zero. The bins next to the central bin (symmetrically to the right and to the left of the central bin) had linear increasing values: 0.01, 0.02, 0.03, and so forth. The maximal value of the weight vector was 1, which corresponded to firing rate values of 2 (in the case of increases) or 0 (in the case of decreases). Firing rates >2 in the PDF were multiplied by 1 in the weight vector to avoid a bias toward increases. Results of the multiplication of the normalized discharge PDF by the weight vector are shown in Fig. 6, *C* and *F*, *right*. To assess the ratio between increases and decreases (the increase/decrease score), we divided the sum of values >1 by the sum of values <1 of the weight-multiplied function.

We repeated the same process with a similar weight vector, which included increasing weight values for firing rates >2 (e.g., in this case a firing rate of 2.01 had a weight of 1.01 instead of 1 as in the usual weight vector). In addition, for both weight vectors we used different smoothing parameters with Gaussian σ values of 20, 50, and 200 ms in addition to 100 ms. We also repeated this analysis by using a firing rate PDF that was normalized to the mean instead of the median.

## APPENDIX C: THE EXTENDED SURPRISE METHOD

Our extended surprise method uses principles similar to those of the surprise burst-detection algorithm (Legendy and Salcman 1985). The method is based on an evaluation of how improbable it is for a certain number of spikes to appear in a time segment of a spike train with a known average firing rate. Our algorithm departs from the original one in two ways: *1*) we use a different definition of the “core segment,” which enables us to detect longer segments of change in firing rate; *2*) we use the same algorithm for detection of decreases in firing rate as well as increases (bursts) in firing rate.

The algorithm begins by binning the spike train into nonoverlapping segments of 100 ms and calculating the mean and the variability (SD) of the firing rate in these segments. We searched for suspected increase/decrease segments (“core segments”) whose firing rate differed by ≥1SD from the mean firing rate of a 100-ms bin. The relatively low threshold of 1SD was chosen to enable a sensitive detection procedure that allowed us later on to test the discharge increases/decrease ratio for a wide range of thresholds. In addition, it prevented bias toward detection of putative increases versus decreases due to the smaller range of firing rate below the mean (the firing rate is bounded from below by zero, whereas the upper bound, given by the refractory period and other cellular parameters, is about 500 spikes/s for GP neurons; Kita and Kitai 1991). The beginning and the end of the core segments were adjusted to the closest spike time before and after the core segment accordingly, both for suspected increases and decreases in firing rate.

After detecting the core segment the program conducts a series of iterative tests in which the current surprise value (reflecting the probability of finding the current number of spikes in the current time duration) is compared with a new surprise value obtained after omitting or adding a single interspike interval to the tested segment. The Poisson surprise value (*S*) of a segment is calculated by (C1) where *P*(*n*) is the probability of finding such a segment assuming a Poisson process. This value is defined for decreases in firing rate as *P _{d}*(

*n*) and for increases in firing rate as

*P*(

_{i}*n*). In cases of decreases in firing rate

*P*(

_{d}*n*) is the probability of finding

*n*spikes or less in an interval of

*T*ms. This probability is based on a Poisson distribution with a mean value of

*rT*, where

*r*is the probability of emitting a spike in a 1-ms bin calculated from the mean firing rate (spikes/s) of the cell and

*T*is the duration of the tested segment. Thus (C2) In cases of increases in firing rate

*P*(

_{i}*n*) is the probability of finding

*n*spikes or more in an interval of

*T*ms. This probability is based on the same Poisson distribution as before. Thus (C3)

This value was calculated in Matlab software by subtracting from 1 the probability of finding *n* − 1 spikes or less in an interval of *T* ms, as follows (C4)

In both cases (increases and decreases) the number of spikes in a segment (*n*) was defined to include the first spike and exclude the last spike in the segment (e.g., for an increase segment with 10 spikes including the beginning spike and the ending spike, *n* = 9; Legendy and Salcman 1985). Note that *Eqs.* C*2* and C*4* include the case where *i* = 0 (zero spikes), although, as described earlier, the minimal number of spikes in an interval is 1. However, including the zero case lowers the surprise value of decreases and raises the surprise value of increases in discharge rate. Therefore as required by our design, it biases the results against decreases in discharge rate.

Matlab's numerical precision of the surprise value was not sufficient in some cases of increases in firing rate (i.e., the cumulative probability of finding the current number of spikes or more was very small under the Poisson assumption and was numerically rounded to zero). In these cases we used the Stirling approximation (C5) In our physiological data, we used this approximation in 349 cases of increases in discharge rate, which equals 0.11% of the total number of increases in firing rate and 0.06% of the total number of changes in firing rate (increases and decreases).

Unlike the original surprise algorithm (Legendy and Salcman 1985), we repeated the maximization process from each side of the core segment both for adding intervals to the core segment as well as for omitting intervals from the core segment. The algorithm begins by finding the number of additional intervals after the core segment which maximizes the surprise value. Then the algorithm calculates the number of omitted intervals from the end of the original core segment which maximize the surprise value. The surprise value of these two conditions (maximization of adding intervals and maximization of omitting intervals) are compared and the condition that has a higher surprise value is chosen. After finding the segment that maximizes the surprise value by adding or omitting intervals at the end of the core segment, we repeated the same process at the beginning of the core segment. Adjacent final segments were merged into one continuous segment only if they were of the same category (increase or decrease in discharge rate) and they overlapped (i.e., the terminating spike of the first segment occurred after the starting spike of the second segment). To avoid bias in detecting increases or decreases, we did not exclude cases in which the same interval took part both in increases and decreases. However, in our physiological database (*n* = 385 cells) these cases were rare: the mean fraction of intervals that took part in both increases and decreases out of the total number of intervals that were attributed to increases or decreases per cell was 0.0003 ± 0.0007.

The time-begin, time-end, and surprise values of all detected increase and decrease segments were stored. We used different surprise values as thresholds for defining a decrease or an increase in firing rate. We first calculated the frequency of increases and the frequency of decreases under different thresholds of surprise values in a single cell. In cases where no increases or decreases were found under a certain surprise threshold, the frequency was defined as zero. We then calculated the ratio between the frequency of increases and the frequency of decreases in each of the surprise thresholds. We ignored surprise thresholds in which the frequency of decreases or increases equaled zero (we did this for decreases to avoid dividing by zero in the calculation of the increase/decrease ratio, and added this exclusion to increases to keep the analysis symmetric). Since the number of segments with increases or decreases in the discharge rate gradually decreases as the surprise threshold rises, the increase/decrease ratio estimate is noisier for higher surprise values. We chose the surprise thresholds as samples out of the actual surprise values of a single neuron for single neuron analysis, or from a data set (i.e., neurons from one of the monkeys recorded from a certain nucleus in a certain state) for population analysis. The sampling process was designed to obtain ≥500 points for each curve uniformly selected out of the full and sorted pool of surprise values of a single neuron or data set. Examples of the increase/decrease ratio curves of single cells are shown in Fig. 9, *B* and *D*. The increase/decrease ratio curves were averaged over all cells in a given data set. This enabled us to obtain a summarized plot of the mean increase/decrease ratio along a range of surprise values in a specific data set (Fig. 11).

## APPENDIX D: SIMULATED NEURONS

We validated the three methods we used (HMM, weighted PDF, and the extended surprise) by applying them to 2,900 simulated neurons of two types: *1*) Poissonian neurons; and *2*) Markov model neurons.

*1*) The neurons in this group were modeled as a realization of a Poisson process with a constant firing probability (*p*) for each time bin (Δ*t*) and an absolute (zero firing probability) refractory period with a length of τ_{ref} bins. The firing probability (*p*) was derived from a normal distribution (μ_{p}, σ_{p}). The simulations were obtained by using the following values: μ_{p} = 0.055, σ_{p} = 0.015, Δ*t* = 1 ms. The values were chosen according to our neural database values (Table 1). We used two values for the refractory period: *1*) τ_{ref} = 0 (no refractory period); *2*) τ_{ref} = 5 ms (absolute refractory period), and in each group we generated 100 cells. The firing probability of cells with an absolute refractory period was adjusted so that the mean firing rate would be similar to the firing rate in the nonrefractory period cells: *P _{ref}* =

*P*/[1 − (

*P*× τ

_{ref})].

These parameters yielded neurons with mean firing rates of 53.07 ± 14.13 spikes/s (in the nonrefractory group) and 53.97 ± 7.68 spikes/s (in the refractory group). The length of the simulated spike trains was 10^{6} bins = 1,000 s.

*2*) The Markov model neurons were generated as a discrete Markov process (Rabiner 1989). The activity of the neuron was assumed to be in one of three states: a baseline state, an increase state and a decrease state. In any time bin (Δ*t* = 1 ms) the neuronal activity was assumed to be in one of several states (1 ≤ *i* ≤ *N*, *N* = 3), with a distinct probability *b _{i}* of emitting a spike in this state. The transition probability from state

*i*at time bin

*t*to state

*j*at time bin

*t*+ 1 is denoted by

*a*. The initial state probability (the probability of being in state

_{ij}*i*in the first time unit) is denoted by π

_{i}. In all simulated neurons we assumed that the neuron was in the baseline state at time unit

*t*= 0; thus π

_{i}=

*A*

_{1i}.

Our simulated database contained neurons with different discharge increase and decrease ratios. In the main balanced condition between discharge increases and decreases, we used the following parameters

These parameters correspond to a three-state model. The mean firing rate in state *i* is derived from *b _{i}* (where

*i*= 1, 2, 3 for background, increase, and decrease states respectively) and equals 50, 80, and 20 Hz at the baseline, increase and decrease states. Thus the mean amplitude change is 30 Hz for both increases and decreases. The mean frequency of increases or decreases is derived from

*a*

_{1i}(where

*i*= 2 and 3 for increase and decrease, respectively) and equals 10/min, in both cases. The average time in the same state (sojourn time) is derived from the probability

*a*of staying in the same state (sojourn probability). Under the Markov-model assumptions, the sojourn times are distributed geometrically, and the average sojourn time in state

_{ii}*i*is 1/(1 −

*a*). In our case this value equals 300 ms for both increases and decreases. Note that we did not allow transitions between increases to decreases or vice versa.

_{ii}Using the Markov model we simulated three groups of “nonbalanced” neurons. In each group we changed a single parameter out of the following: the mean frequency of discharge increases and decreases, the mean duration of discharge increases and decreases, the mean firing rate change (amplitude change) of discharge increases and decreases. In each group we used three possible values for discharge increases and decreases; thus the number of combinations in a group equaled 9. In the first group (frequency) we used the following parameters: 5, 10, and 20/min. In the second group (duration) we used the following parameters: 100, 300, and 500 ms. In the third group (amplitude change from baseline) we used the following parameters: 10, 30, and 50 Hz. For example, in the first group the possible combinations were: 5–5, 5–10, 5–20, 10–5, 10–10, 10–20, 20–5, 20–10, 20–20 (discharge increases–decreases/min). These simulation values were chosen to cover the frequency range of spontaneous pauses of GPe neurons (Elias et al. 2007) and the typical response amplitude and duration of HFD pallidal and SNr neurons (Jaeger et al. 1995; Sato and Hikosaka 2002; Turner and Anderson 1997). As a convention in this study, the first number represents the increase parameter and the second number represents the decrease parameter. The 27 combinations described earlier do not include interactions between two or three parameters (e.g., changes of both frequency and duration). However, this array enables a basic evaluation of different methods without excessive complexity.

Using the Matlab random number generator we generated 100 cells with the same simulation parameters. Thus our total Markovian database consisted of 3 groups × 9 configurations/group × 100 cells/configuration = 2,700 simulated cells. The length of each simulated neuron was 10^{6} bins = 1,000 s.

## GRANTS

This study was partly supported by the Hebrew University Netherlands Association (HUNA) “Fighting against Parkinson” grant to H. Bergman, a Foulkes Foundation fellowship to S. Elias, and a fellowship from Myers-JDC-Brookdale Institute of Gerontology and Human Development and ESHEL, the Association for the Planning and Development of Services for the Aged in Israel to S. Elias.

## Acknowledgments

We thank O. Marmor for assisting in recording, G. Heimer and J. A. Goldberg for sharing their data, G. Goelman for MRI, V. Sharkansky and S. Freeman for technical assistance, and E. Singer for language editing.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society