Thalamocortical Relay Fidelity Varies Across Subthalamic Nucleus Deep Brain Stimulation Protocols in a Data-Driven Computational Model

Yixin Guo, Jonathan E. Rubin, Cameron C. McIntyre, Jerrold L. Vitek, David Terman


The therapeutic effectiveness of deep brain stimulation (DBS) of the subthalamic nucleus (STN) may arise through its effects on inhibitory basal ganglia outputs, including those from the internal segment of the globus pallidus (GPi). Changes in GPi activity will impact its thalamic targets, representing a possible pathway for STN-DBS to modulate basal ganglia-thalamocortical processing. To study the effect of STN-DBS on thalamic activity, we examined thalamocortical (TC) relay cell responses to an excitatory input train under a variety of inhibitory signals, using a computational model. The inhibitory signals were obtained from single-unit GPi recordings from normal monkeys and from monkeys rendered parkinsonian through arterial 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine injection. The parkinsonian GPi data were collected in the absence of STN-DBS, under sub-therapeutic STN-DBS, and under therapeutic STN-DBS. Our simulations show that inhibition from parkinsonian GPi activity recorded without DBS-compromised TC relay of excitatory inputs compared with the normal case, whereas TC relay fidelity improved significantly under inhibition from therapeutic, but not sub-therapeutic, STN-DBS GPi activity. In a heterogeneous model TC cell population, response failures to the same input occurred across multiple TC cells significantly more often without DBS than in the therapeutic DBS case and in the normal case. Inhibitory signals preceding successful TC relay were relatively constant, whereas those before failures changed more rapidly. Computationally generated inhibitory inputs yielded similar effects on TC relay. These results support the hypothesis that STN-DBS alters parkinsonian GPi activity in a way that may improve TC relay fidelity.


The delivery of high-frequency stimulation to the subthalamic nucleus (STN) or other target areas, through a surgically implanted electrode, has become a widely used therapeutic option for the treatment of Parkinson's disease (PD) and other neurological disorders (Benabid et al. 2006). The mechanisms underlying the effectiveness of deep brain stimulation (DBS), however, remain unclear and under debate. Multiple studies have shown that pathological rhythmicity emerges in certain subsets of cells within the basal ganglia in parkinsonism (Bergman et al. 1994; Brown et al. 2001; Hurtado et al. 1999, 2005; Levy et al. 2003; Magnin et al. 2000; Nini et al. 1995; Raz et al. 2000). Therefore DBS for PD may work by eliminating or modifying such pathological signals. Initial attempts to address this concept focused on the possibility that DBS blocks neural activity, creating a physiologic lesion (Beurrier et al. 2001; Filali et al. 2004; Magarinos-Ascone et al. 2002; Tai et al. 2003; Welter et al. 2004). According to this theory, suppression of thalamic firing by inhibition from basal ganglia output areas, such as the pallidum, is reduced by DBS, and through this reduction DBS restores the capability of the thalamus to engage in appropriate movement-related activity (Benabid et al. 2001; Benazzouz et al. 2000; Obeso et al. 2000; Olanow and Brin 2001; Olanow et al. 2000).

Recent experimental and computational results, however, suggest that neurons directly downstream from stimulated regions may in fact be activated by DBS (Anderson et al. 2003; Hashimoto et al. 2003; Hershey et al. 2003; Jech et al. 2001; McIntyre et al. 2004; Miocinovic et al. 2006; Paul et al. 2000; Windels et al. 2000, 2003). These results support the alternative idea that DBS works by replacing pathological rhythms with regularized firing activity (Foffani and Priori 2006; Foffani et al. 2003; Garcia et al. 2005; Grill et al. 2004; Meissner et al. 2005; Montgomery and Baker 2000; Vitek 2002). In past theoretical work, we offered a computational implementation of this idea (Rubin and Terman 2004). We used Hodgkin-Huxley-type models of cells in the indirect pathway of the basal ganglia (Terman et al. 2002) to generate inhibitory output trains, which served as synaptic inputs to a model thalamocortical (TC) relay cell. In this previous model system, we assessed TC cell activity under stereotyped representations of normal, parkinsonian, and DBS conditions. Our simulations and analysis demonstrated and explained a mechanism by which pathological oscillatory or bursty inhibition from the internal segment of the globus pallidus (GPi) to TC cells could compromise the fidelity of TC relay of excitatory signals, whereas elimination of the oscillations within this inhibition, even at levels that are elevated relative to normal conditions, could restore TC cells’ relay capabilities (Rubin and Terman 2004).

In this study, we use GPi spike trains recorded from normal control monkeys and from parkinsonian monkeys (Hashimoto et al. 2003), with or without DBS of the STN region, as the source of inhibitory inputs to our model TC cells. By doing so, we circumvent the controversy surrounding the effects of DBS at the stimulation site. Within this theoretical framework, we are able to test how biologically observed changes in GPi neuronal activity affect TC signal transmission, both in a single-model TC cell and in a heterogenous population of model TC cells. TC relay fidelity is evaluated using a train of external excitatory stimuli applied to the same model TC cells that receive the recorded inhibitory synaptic inputs from GPi. Our results show that there is a significant decline in the ability of the TC cells to relay the excitatory stimuli when they are exposed to GPi signals recorded under parkinsonian conditions in the absence of DBS or with sub-therapeutic DBS, defined by its failure to induce a therapeutic effect on motor symptoms, relative to GPi data recorded from normal monkeys. Moreover, relay effectiveness is restored to nonparkinsonian levels by GPi signals recorded under parkinsonian conditions in the presence of therapeutic DBS, which induced a measurable improvement in motor symptoms. Interestingly, while response failures across a population of TC cells tend to occur on similar trials in the parkinsonian and sub-therapeutic cases, failures occur asynchronously under therapeutic STN DBS as well as under normal conditions, which would moderate their downstream effect. Finally, to extend these results, we harness a purely computational approach that allows us to systematically vary the rhythmicity and degree of correlation within the inhibitory inputs that TC cells receive. Our results show that moderately increasing the burstiness and correlation of inhibitory spike trains, as might be expected in a transition from normal to parkinsonian conditions, leads to a gradual loss of relay fidelity, while a further transition to tonic high-frequency, highly correlated inhibitory signals, as may occur in clinically effective DBS (Hashimoto et al. 2003), leads to significant restoration of effective relay.


Proposed mechanism for DBS effectiveness

In awake states, TC cells serve to relay excitatory inputs (Steriade et al. 1997). The TC population targeted by GPi cells likely is involved in the relay of excitatory inputs between cortical areas (Guillery and Sherman 2002a,b; Haber 2003). The basic idea being explored in this paper is that changes in inhibitory output from the GPi to its target TC cell population affect the relay reliability of these TC cells, defined in terms of the generation of TC activity patterns that match the inputs to TC cells. Specifically, parkinsonian conditions induce oscillations, burstiness, and enhanced correlations in GPi outputs, and these effects are hypothesized to compromise relay fidelity. We further hypothesize that the effectiveness of DBS is due to the replacement of pathological GPi firing patterns with more regular activity. While this regular activity may in fact be overly regular, and may occur at a higher frequency, relative to the activity that occurs in nonparkinsonian states, it nonetheless restores thalamocortical relay reliability. This concept is illustrated schematically in Fig. 1.

FIG. 1.

Hypothesized mechanism for deep brain stimulation (DBS) effectiveness. In each of the 3 cases shown, the target thalamocortical (TC) cell receives inhibitory inputs from the internal segment of the globus pallidus (GPi), which affects its relay of an excitatory drive. In the normal case, the inhibition is irregular and relatively weak due to low correlation levels (represented by ░), and the TC cell successfully relays its inputs. In the parkinsonian case, inhibition is more bursty and stronger (▪) due to enhanced correlations. During each inhibitory burst, the TC cell fails to respond to its drive (i.e., misses), while its response is excessive (i.e., bad) between bursts. In the case of supraclinical or therapeutic DBS, inhibition is strong but quite regular. Despite the strength of the inhibition, successful TC relay is restored.

These effects on TC relay in parkinsonian and DBS conditions remain to be demonstrated experimentally, but they were shown to arise in a previous, purely computational study (Rubin and Terman 2004) where a possible dynamical mechanism that could yield these results was also explained. The fundamental hypothesis from our original study was that DBS leads to tonic, regular inhibitory input to the TC cells, and this allows the activation and inactivation levels of TC cell membrane ionic currents to equilibrate, such that reliable relay can occur, as long as excitatory inputs are not excessively rapid. During parkinsonian conditions, the inhibitory output of GPi features synchronized oscillations with bursting activity. When a significant increase in the level of inhibition of TC cells associated with such oscillations occurs, a period of re-equilibration of the TC ionic currents ensues. During this time, it is difficult for the TC cells to reliably respond to excitation (Jahnsen and Llinas 1984a). Further, after currents have equilibrated to a high level of inhibition, a relatively abrupt decrease in inhibition can lead to an excessive or bursty TC response to excitation due to increased availability of spike-generating and -sustaining currents (Jahnsen and Llinas 1984a,b). We propose that therapeutic DBS reduces this oscillatory activity in GPi and TC cells, thereby improving the ability of TC cells to relay information.

Model TC cells

The model used for the TC cells is a slightly modified version of that used in our earlier study (Rubin and Terman 2004), which is itself a simplification of an earlier formulation (Sohal and Huguenard 2002; Sohal et al. 2000). In this model, the current-balance and ionic activation equations take the form Math Math Math(1)

In the preceding equations, the terms IL = gL[vEL], INa = gNam3(v)h[vENa], and IK = gK[0.75(1 − h)]4[vEK] are leak, sodium, and potassium spiking currents, respectively, with square brackets denoting multiplication. Note that we use a standard reduction in our expression for the potassium current, which decreases the dimensionality of the model by one variable (Rinzel 1985). The current IT = gTp2(v)r[vET] is a low-threshold calcium current. For these intrinsic currents, the forms of the functions and the values of the parameters used appear in Table 1. Note that reversal potentials are given in mV, conductances in mS/cm2, and time constants in ms. Further, we have scaled the parameters such that the capacitance is Cm = 1 μF/cm2. Finally, the resting potential, spike threshold, and responsiveness of the model TC cell, in the absence of inputs, are robust to changes of ionic conductances in the model. Durations of rebound bursts, after release from hyperpolarizing input, may jump abruptly by tens of milliseconds as gT is varied, however, when an additional spike is appended to the burst. As is typical for conductance-based models, the model is less robust to changes in the threshold and slope constants within its nonlinear terms; however, its robustness is comparable to other models of this type.

View this table:

TC cell model functions and parameters

Additional terms in Eq. 1 refer to inputs to the TC cell model. The equations and parameter values relevant to these terms are summarized in Table 2 with the same units used as in Table 1. Iext corresponds to a constant background input, chosen at Iext = 0.44 nA/cm2 to yield a firing rate of roughly 12 Hz in the absence of other inputs and held fixed at this level throughout all simulations. The value chosen places the model TC cell near transition from silent to spontaneously oscillatory in the absence of synaptic inputs. Similar results were obtained whether the model TC cell was intrinsically silent or oscillatory. By choosing Iext near the transition point, we achieved wide variations in TC cell behaviors when we introduced variability into the set of model TC cell parameters as discussed in the following text. IGi→Th denotes the inhibitory input current from the GPi and will also be discussed in the following text. IE represents simulated excitatory synaptic signals to the TC cell. We assume that these are sufficiently strong to induce a spike (in the absence of inhibition) and therefore may represent synchronized inputs from multiple presynaptic cells. In the model, IE takes the form gEs[vvE] where gE = 0.05 μS/cm2, so that maximal input is super-threshold, where vE = 0 mV, and where Math(2)

View this table:

Inputs to the TC cell

In Eq. 2, α = 0.8 ms−1 and β = 0.25 ms−1. Because we do not have an explicit representation of a presynaptic neuron in the model, we use the function exc(t) to control whether the excitatory input is on or off. Specifically, exc(t) = 1 during each excitatory input, whereas exc(t) = 0 between excitatory inputs. We used two general forms of time course for the binary signal exc(t), namely periodic and stochastic, as done in past work (Rubin and Terman 2004). In the periodic case Math where the period p = 50 ms and duration d = 5 ms, and where H(x) is the Heaviside step function, such that H(x) = 0 if x < 0 and H(x) = 1 if x > 0. That is, exc(t) = 1 from time 0 up to time d, from time p up to time p + d, from time 2p up to time 2p+d, and so on. A baseline input frequency of 20 Hz is consistent with the high-pass filtering of corticothalamic inputs observed in vivo (Castro-Alamancos and Calcagnotto 2001); at this input rate, the model TC cells rarely recover and fire spontaneous spikes between inputs, which simplifies our analysis. In the stochastic case, input times are selected from a Poisson distribution, with an enforced pause of 20 ms between inputs to avoid excessive firing, with the same input duration and amplitude as in the periodic case and with a mean input frequency of 20 Hz. In simulations done with stochastic inputs, results shown represent averages over five simulations, each with a different random input pattern. The use of a stochastic excitatory input provides one measure of the robustness of our results to noise. In some simulations, a small-amplitude white-noise term ξ is also included in the voltage equation.

The choice of gE was motivated by the conjecture that strong inputs would represent important signals and that differences in TC relay of strong inputs would have the most significant impact on downstream processing. At the same time, it is unlikely that even strong inputs would be perfectly synchronized. The values of the rate parameters α, β, and d were selected based on corticothalamic excitatory inputs recorded in vivo (Castro-Alamancos and Calcagnotto 2001), assuming that IE represents a set of temporally proximal, but imperfectly aligned, cortical inputs to a TC cell. Our qualitative results are robust to variations in these parameters.

In one set of simulations, we feed the same input currents IGi→Th and IE into all members of a heterogeneous population of model TC cells. In the absence of experimental data on the variability of particular conductances within TC cells, we chose to form the heterogeneous population by selecting gNa, gL, and gT from normal distributions with means given by the values in Table 1 and with SDs given by 20% of these values, which were sufficiently large to yield a wide variation in intrinsic TC spike frequencies in the absence of inputs, without major changes in most other spike-related characteristics, as discussed in the text following Eq. 1.

Experimentally obtained GPi data

Single-unit extracellular recordings of neurophysiologically identified GPi neurons were acquired with glass-coated platinum-iridium microelectrodes in three rhesus macaques (Macaca mulatta). One animal was a normal nonparkinsonian control, and two animals were rendered parkinsonian with 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) via a single injection through the internal carotid artery (Hashimoto et al. 2003). The parkinsonian animals developed a stable disease state characterized by contralateral rigidity and bradykinesia and had a chronic DBS electrode implanted in the STN region (Hashimoto et al. 2003). The chronic stimulating electrode was connected to a programmable pulse generator (Itrel II, Medtronic) implanted subcutaneously in the monkey's back. The stimulating lead was a scaled-down version of the chronic stimulation electrode used in humans (Model 3387, Medtronic). The cylindrical lead consisted of four metal contacts each with a diameter of 0.75 mm, height of 0.50 mm, and separation between contacts of 0.50 mm. The most effective pair of electrode contacts in the STN region was chosen for bipolar stimulation in each animal after evaluation of the clinical effects of the stimulation (Hashimoto et al. 2003). In both the normal and parkinsonian monkeys, spontaneous neuronal activity (with the animal at rest and the head fixed) of electrophysiologically identified GPi neurons was recorded. In the parkinsonian monkeys, GPi activity was also recorded during DBS of the STN region. Stimulation parameters were selected to address two conditions in each animal: stimulation parameters that produced therapeutic benefit and stimulation parameters subthreshold for a therapeutic effect. The therapeutic effectiveness of DBS was assessed with two quantitative measures of bradykinesia as well as a subjective evaluation of rigidity provided by a trained rater. In each animal, therapeutic stimulation settings were determined, and then sub-therapeutic settings were obtained by reducing stimulus amplitude until therapeutic benefit was no longer detected (Hashimoto et al. 2003). Specifically, DBS was applied at a frequency 136 Hz with therapeutic benefit obtained at 3.3 or 1.8 V and a pulse width of 90 or 210 μs, depending on the animal, and subthreshold stimulation at 2 or 1.4 V, again depending on the animal. To analyze neural activity during stimulation, a template of the stimulus artifact was constructed by averaging across all peristimulus segments. The stimulus artifact template was then subtracted from the individual traces, and neuronal spikes were detected (Hashimoto et al. 2002, 2003).

Collections of several cells from each of the three animals were used in the analysis. The cells were selected from a database of recordings to be representative of the population in the normal, parkinsonian, sub-therapeutic DBS, and therapeutic DBS cases. Three general characteristics were used to select the cells. First, the experimental recording had good to excellent isolation of the single unit. Second, the average firing rate of the unit closely corresponded to the average population firing rates for GPi cells in the four respective cases (Hashimoto et al. 2003; Wichmann et al. 2002). Finally, the coefficient of variation of the firing rate was used to identify cells with firing patterns representative of the four respective cases. The particular firing characteristics of the cells used are summarized in Table 3 with relevant values from the literature provided for comparison.

View this table:

Firing characteristics of GPi cells

Inhibitory inputs to TC cells, derived from GPi data

In most simulations, we used experimentally recorded data, as discussed in the preceding text, to represent the GPi spike times. For systematic exploration of the effects of particular features of the inhibitory input, however, we used computationally generated GPi spike times. In both cases, the synaptic inhibition from the GPi to a single model TC cell in our simulations took the form Math(3) where the summation is over the synaptic activation variables sj of the presynaptic GPi cells, and where the inhibitory synaptic reversal potential Esyn = −85 mV (Lavin and Grace 1994) and synaptic conductance gsyn = 0.066 μS/cm2. At each spike time of the corresponding GPi cell, the variable sj was reset to 1, after which it decayed via the equation Math(4) with βinh = 0.04 ms−1. We used a relatively large synaptic conductance and a synaptic decay rate that is somewhat slower than that typically found for GABAA-mediated synaptic transmission to make our single input train more representative of multiple, imperfectly synchronized synaptic inputs; this approximation will be improved in future work as multi-unit GPi data are collected experimentally.

Experimental GPi data were recorded from parkinsonian monkeys before, during, and after the application of DBS (Hashimoto et al. 2003). When we used non-DBS and DBS recordings from the same cell, we only used non-DBS recordings from the period before the application of DBS, not from the period after the cessation of DBS, to avoid any residual effects of DBS on GPi neuronal activity. Moreover, we selected data segments by counting back in time from the end of the DBS period, always stopping ≥2 s away from the start of DBS to minimize the possibility of our results being affected by transients associated with DBS onset.

Error index: a measure of TC relay fidelity

The computations in this paper were performed using customized codes simulated in XPPAUT (Ermentrout 2002; see∼phase) and Matlab (The MathWorks, Natick, MA).

We compute an error index to measure the fidelity with which the TC cells respond to excitatory inputs, similar to the error index described previously (Rubin and Terman 2004). Note that we use relay fidelity to refer to the faithfulness of relay such that a TC cell that generates a spike train that is very similar to its input train has achieved high degree of relay fidelity. We are not using fidelity to refer to the generation of similar responses to multiple presentations of the same stimulus, which is a form of reliability considered in some other studies. In brief, the error index that we use consists of the total number of errors divided by the total number of excitatory inputs. Errors can take the form of bad responses or misses. Specifically, for each excitatory stimulus, we record a miss if no TC spike occurs within a designated detection time window after the input arrives. If more than one TC spike occurs during this window, then we record a bad response. Finally, if exactly one TC spike occurs during the window, then we record a bad response if there are one or more additional TC spikes before the next input arrives (see Fig. 1). This algorithm counts at most one error per input; for example, if a TC cell fires multiple spikes after a single excitatory input, then this is just counted as a single bad response. In summary, the error index is given by Math(5) where b denotes the number of excitatory inputs leading to bad responses, m the number of excitatory inputs leading to misses, and n the total number of excitatory inputs. We use a detection window of 10 ms to allow for delays from threshold crossing to action potential generation. Thus an error index of 0 results if one TC spike occurs within 10 ms of each excitatory input and no other TC spikes occur until the next input, corresponding to optimal relay fidelity. With a shorter detection window, some formerly successful responses would be classified as misses. However, we did not observe any bias toward shorter or longer response latencies in any particular inhibition regime, and indeed, we obtained qualitatively similar results in simulations with detection windows of 6 and 12 ms.

In theory, our error index could be susceptible to “false positives,” in which single spike TC responses occur close in time to excitatory inputs, but not caused by the excitatory inputs. Thus as mentioned earlier, we use excitatory input rates that are sufficiently high such that in normal conditions, TC cells rarely recover and fire spontaneous spikes between inputs. Finally, note that our error index gives a direct and straightforward measure of relay success that is well suited for our computational experiments, in which the simplicity of our simulated excitatory inputs and of the relay process does not warrant analysis with more standard, yet more complex and indirect, information theoretic measures.

Burstiness and correlation of inhibitory GPi signals

Much of our analysis concerns ways in which the error index depends on the burstiness and correlation of the inhibitory GPi signals sj. To quantify burstiness, we first perform a simple detection algorithm for high-frequency spiking episodes (HFE). In this approach, we detect all spikes that are preceded by a silent period of ≥12 ms. Each such spike is considered to be the start of an HFE if the next spike follows it by <8 ms. Each subsequent spike is counted as part of the HFE if and only if it occurs within 8 ms of its predecessor. The duration of the HFE is the time from the first spike in the HFE to the last spike in the HFE; see Fig. 2D. More involved statistical methods exist to compensate for chance epochs of high-frequency spikes that fit within a given set of HFE criteria of the type given here (Legendy and Salcman 1985); however, because all such HFE generate similar inhibitory inputs to TC cells in our simulations, there is no reason to try to classify them for the purposes of our study. From the HFE, we compute the elevated spiking time (EST), which is simply the fraction of the simulation time during which HFE occur. Hence, when the EST is zero, the GPi signal consists of low-frequency isolated spikes, a moderate EST corresponds to a highly bursty signal, and a signal with a higher EST is dominated by HFE, corresponding to relatively tonic high-frequency firing.

FIG. 2.

Examples of inputs from GPi cells to TC cells. AC: examples of the high-frequency burst portions from computationally generated GPi signals. A: a low burst rate rb and no overlaps (shared wij) were used to generate these signals, and correspondingly, there are relatively few bursts, leading to a mean elevated spiking time (EST) of 0.14 across the 2 cells. Moreover, the amount of time during which the traces simultaneously exhibit high-frequency firing is small, yielding a small correlation time of 0.082. B: a moderate burst rate rb and 2 overlaps were used to generate these signals, and correspondingly, each GPi trace shows high-frequency oscillations for about half of the total simulation time, with a mean EST of 0.61. Although the times at which these occur are somewhat correlated, due to the overlaps and chance, the fraction of the total simulation time during which the traces simultaneously exhibit high-frequency firing is <1/2 with a correlation time of 0.45. C: with a high rb and 2 overlaps, each trace exhibits high-frequency oscillations for most of the simulation time, yielding an EST of 0.85, and the fraction of time during which the traces simultaneously show high-frequency activity is much closer to 1, yielding a correlation time of 0.73. D: illustration of the algorithm for detection of coincident high-frequency episodes (HFE), applied to experimental data. The times at which HFE occur are read off of GPi spike trains (top 2 panels; HFE times are indicated with thick black segments in all panels). Next, HFE times are compared and times when both cells are engaged in HFE are captured (bottom panel).

The aspect of the correlation between pairs of GPi signals that is most relevant for our study is the temporal relationship of the HFE across the two signals. To obtain a single number that represents this relationship over a simulation of duration T ms, we simply sum the durations of all epochs during which both GPi cells are engaged in high-frequency spiking simultaneously and divide by T, yielding a number between 0 and 1 (Fig. 2D).

Event-triggered averaging, sorted by TC cell responses

An additional computational procedure that we performed on GPi data was event-triggered averaging. In this procedure, we classified excitatory inputs into those that were immediately followed by a missed, a bad, or a successful (i.e., neither missed nor bad) TC response. For each excitatory input that led to a missed response, we extracted a 25-ms segment of the GPi input signal to the TC cell, extending from 20 ms before the start of the excitatory input to 5 ms after its start. The entire time course of each signal was normalized by subtracting off the signal's initial amplitude. We summed these normalized, “miss-triggered” GPi signals and divided by the number of missed responses to generate a miss-triggered average GPi signal. Next we repeated this procedure for bad responses and successful responses to generate a bad-triggered average GPi signal and a success-triggered average GPi signal, respectively. In this averaging process, we combined inhibitory signals leading to the same type of response from all four inhibitory input regimes (nonparkinsonian, parkinsonian without DBS, parkinsonian with sub-therapeutic DBS, and parkinsonian with therapeutic DBS), after verifying that similar signals emerged in all cases. In total, 40 blocks of GPi data, each of 5-s duration, were used. These yielded 280 bad responses, 667 missed responses, and 2,223 successful responses, all of which were included in the averages computed.

Plots of average GPi signals do not include error bars. We chose to omit them because the error bars for averages of GPi signals could be large, despite a very high degree of qualitative similarity, such as when each signal showed an abrupt increase at some time within a given time window, but the precise increase times were rather diverse. A similar issue arises in averaging the action potential responses of a neuron over multiple stimulus presentations or in multi-unit recordings, in averaging over action potentials generated by different cells in response to the same stimulus (e.g., Kapfer et al. 2007). Following the procedure used by Kapfer et al. (2007), instead of plotting error bars, we complement plots of average signals with data from a sample of individual signals that contributed to the averages, selected completely at random.

Jittered inputs

Note that the experimental GPi data used in this study consist of single-unit recordings acquired with a single electrode. Therefore it was not possible to use this data directly to explore how correlations among multiple GPi inputs to TC cells contribute to the TC cell relay fidelity. Because we did not have this option, in some simulations, we used the single-unit GPi recordings to generate multiple GPi signals to each TC cell. To do this, we first formed N identical copies of a single GPi spike train. We indexed the spike times within this train as {t1, t2,…, tp}. Next, we introduced jitter by selecting values ξij, i = 1,…, N, j = 1,…, p, from a normal distribution with amplitude σ (see results for σ values used). These were used to form the new spike trains {t11, t12,…, t1p},…,{tN1, tN2,…, tNp} with tij = tj + ξij. After some experimentation, we found that the qualitative trends induced by this jittering process are already apparent with N = 2. Given this observation, we restrict our results to the case of N = 2, and we also turn to simulated GPi inputs to explore more thoroughly the effects of different activity patterns and different levels of correlations among inhibitory signals.

Computational GPi inputs and their burstiness and correlation

By using purely computational GPi input signals, we were able to explore systematically how changes in input ESTs and the degree of correlation between inputs affect TC relay. For simulated GPi inputs, each signal sj, j = 1,2, in Eq. 2 was formed using a computational procedure, rather than using experimental data, based on a combination of five independent point processes, wij, i = 1,…, 5; see Fig. 7. Each point process wij was produced by a set of four Poisson processes. One Poisson process (p1) was used to generate isolated spike times. A set of three additional Poisson processes were used to generate bursts of high-frequency activity that were superimposed on the isolated spikes. Specifically, a primary Poisson process (p2) selected HFE onset times with rate rb, while within bursting HFEs, a secondary process (p3) produced spike times, at high frequencies. Finally, HFE durations were selected randomly from a third, independent Poisson process (p4), with a minimum duration of 10 ms and a mean duration of 25 ms, for all rb. For each GPi cell in the computational case, the EST was computed as the sum of the durations of all HFEs for the point processes used to form the signal sj for that cell. This approach is computationally simpler than basing the EST on particular spike times and interspike intervals within each HFE, as was done in the experimental case, although it yields EST values that are systematically larger than those obtained in the experimental case.

The five point processes wij were used to generate a single continuous time input signal sj(t) (see Figs. 2 and 3). Specifically, at each spike time within any of the wij, the variable sj(t) was reset to 1, after which it decayed continuously via Eq. 4. This approach, of generating a continuous time signal sj(t) from a collection of point processes wij, allows for parametric control of the degree of burstiness and the spike rate of each wij, and hence of each sj(t) (Tateno and Robinson 2006). The reason that we used multiple signals wij for each sj(t) is that this allowed us to control the correlation across the sj by using some of the same signals wij for different j (Galan et al. 2006; see Fig. 7). We refer to the number of signals wij shared by two GPi cells as the number of overlaps between them.

FIG. 3.

TC relay fidelity improves with clinically effective DBS. AD: the central trace in each plot shows voltage vs. time for the model TC cell. The voltage scale on each plot applies to this trace. Offset above each such trace, experimentally recorded GPi spike times (discrete events) are shown along with the inhibitory signal s1 that these spike times are used to generate (continuous curve, above the spike times, shows s1, with amplitude scaled 100-fold for visibility). Offset below each TC voltage trace, simulated excitatory input signals are shown (scaled by a factor of 3 for visibility). Note that the same excitatory input signals were used for all examples shown here and that TC spikes may lag excitatory input times by a few milliseconds, corresponding to delays from threshold crossing to spike generation. A: control (nonparkinsonian); EST = 0.05. B: parkinsonian without DBS; EST = 0.15. C: parkinsonian with sub-therapeutic DBS; EST = 0.27. D: parkinsonian with therapeutic DBS; EST = 0.55. E and F: error index against EST calculated from simulations of 5-s blocks of data from all 4 cases. In these plots, results for the different cases are color coded (purple: normal, 2 blocks from each of 3 cells; blue: parkinsonian without DBS, 3 blocks from each of 3 cells and 4 blocks from 1 cell; green: parkinsonian with sub-therapeutic DBS, 6 blocks from 1 cell and 2 blocks from 1 cell; red: parkinsonian with therapeutic DBS, 6 blocks from 1 cell and 5 blocks from another cell). Across the 3 parkinsonian cases, each symbol corresponds to the use of data from a particular GPi cell. For example, results indicated by a blue diamond and a red diamond were obtained using data from the same GPi cell, recorded in the absence of DBS and with therapeutic DBS, respectively. E: results from 20-Hz periodic excitatory inputs. F: results from excitatory inputs generated by a Poisson process with a minimum time interval of 20 ms imposed between inputs.

To form the total synaptic input conductance to the TC cell as a function of time, the signals s1(t), s2(t) were summed, as indicated in Eq. 3, and multiplied by gsyn = 0.04 μS/cm2. This maximal synaptic conductance value is smaller than was used in the experimental case to compensate for the replacement of a single experimental GPi signal with a pair of computational ones.


With experimentally obtained GPi inputs, clinically effective DBS improves TC relay fidelity

We generated GPi inputs to our model TC cell using GPi spike trains obtained from experimental recordings from a normal monkey as well as from two parkinsonian monkeys in the absence of DBS, during sub-therapeutic DBS, and during therapeutic DBS (Hashimoto et al. 2003), as described in methods. In each simulation, a single GPi train was used, and hence the sum ∑j sj in Eq. 3, became simply s1. Figure 3 (top traces in each panel) shows typical examples of the experimentally recorded GPi spike times and the resulting GPi signal s1 from each regime. The pattern of GPi activity recorded in parkinsonian conditions in the absence of DBS led to a GPi signal (Fig. 3B, top trace) that is much more phasic, featuring jumps between high and low states, than the relatively constant signal that appeared when therapeutic DBS was present (Fig. 3D, top trace) or in nonparkinsonian conditions (Fig. 3A, top trace). In each case, an excitatory input train was delivered (Fig. 3, bottom traces) and the effectiveness of the TC cell at relaying this train was assessed.

If perfect relay fidelity were achieved, the TC cell would exhibit one voltage spike for each input pulse, possibly with a short lag due to the delay between threshold crossing and actual spiking. In the absence of DBS and with subclinical DBS, however, the TC cell failed to respond to many of the inputs and generated bursts of multiple spikes to other inputs (Fig. 3, B and C, middle traces). These results contrast strongly with the normal and therapeutic DBS cases, which show a near-perfect relay performance (Fig. 3, A and D, middle traces).

We calculated the error index based on the computational TC cell's relay performance for each of the four cases, namely control, PD (no DBS), sub-therapeutic DBS, and therapeutic DBS, both for periodic excitation and for stochastically timed excitation. Results are shown in Fig. 3, E (periodic) and F (stochastic), where each data point represents 5 s of simulation time, with nonoverlapping 5-s GPi data segments used, and is plotted as a function of the EST of the GPi input signal, computed as described in methods. It is important to note that for the GPi recordings involved, all available data were used; that is, we did not select out particular simulation periods based the resulting error indices. The values of the error index show that TC cell relay success depends strongly on which form of inhibitory input the cell receives. Indeed, in both the periodic and the stochastic excitation cases, the mean performances across the four regimes were statistically significantly different (ANOVA, P < 0.0001), and a posteriori pairwise comparisons yielded significant differences across all pairs of regimes in both cases as well (Tukey's honestly significant difference, P < 0.01 for all pairs), except no significant difference was found between therapeutic DBS and normal cases either with periodic excitation or with stochastic excitation. Similar results were obtained with variations in the rise and decay times of the excitatory input signals and in the detection window used to define successful TC responses as well as with the introduction of small noise as shown in Eq. 1. Once rise times dropped by ∼20%, the statistical significance of the differences in error index values between some cases, particularly sub-therapeutic DBS/normal, began to degrade; however, the qualitative distinctions between these values remained.

Differences in GPi signals precede different TC cell responses

As described in methods, the TC cell response to each excitatory input was classified as a miss if the TC cell failed to spike within a prescribed time window following the input, a bad response if the TC cell generated multiple spikes in response to the input, or a successful response. Misses and bad responses raised the error index, while successful responses did not. All three types of responses were found, in differing proportions, in the four scenarios of normal and of parkinsonian with no DBS, with sub-therapeutic DBS, and with therapeutic DBS. To analyze further the way in which the inhibitory signal to the TC cell contributed to its responses, we performed the averaging procedure described in methods on the same GPi input signals used to compute the error index scores (Fig. 3, E and F). We observed important differences across the resulting miss-, bad-, and success-triggered GPi signals (Fig. 4; n = 667 miss, n = 220 bad, n = 2223 success). GPi inputs that preceded TC cell misses showed a substantial rise in strength over the 25-ms time interval considered. In the face of such a rise in inhibition, the TC cell would require additional deinactivation of its spike-generating currents, namely INa and IT in the TC model (1), relative to their resting levels, to respond to an incoming excitatory stimulus (Rubin and Josic 2007; Rubin and Terman 2004). This deinactivation occurs relatively slowly, however, and thus would typically require more than the 25 ms available here.

FIG. 4.

Average GPi signals preceding different types of TC cell responses to excitatory inputs are qualitatively different. A: the 3 traces shown were formed by averaging over 25-ms segments of normalized GPi signals s1, spanning the arrival times of excitatory inputs to a TC cell. The signals were aligned such that the excitatory input arrival times occurred at 20 ms as indicated by the vertical dashed line in the figure. Signals were averaged separately for excitatory inputs that produced TC cell misses (n = 667), bad responses (n = 280), or successful responses (n = 2223). B–D: the values at 0 ms, 20 ms (i.e., excitatory input arrival time), and 25 ms for a randomly selected sample of 10 normalized miss-triggered (B), bad-triggered (C), or success-triggered (D) signals, from the sets of signals used to generate the averages shown in A.

Conversely, GPi inputs that preceded bad TC cell responses showed a substantial decline in strength over the 25-ms interval considered. Recall that what we classify as bad responses consist of multiple spikes fired in response to single excitatory inputs because such responses do not reflect the content of the input signals. In the presence of a strong inhibitory input, deinactivation of a TC cell's spike-generating currents will occur. The resulting enhanced availability of these currents will allow for successful responses in the presence of sustained inhibition. When followed by a relatively rapid drop in inhibition, however, as seen in the bad-averaged signal in Fig. 4, this additional deinactivation will lead to an excessive response to excitatory inputs (Rubin and Terman 2004) until it can be negated by a subsequent slow inactivation of the currents involved.

Finally, GPi inputs that preceded successful TC cell responses were relatively constant and therefore avoided the generation of current imbalances. Interestingly, the roughly constant averaged inhibition level in this case was relatively high (data not shown). This is consistent with the notion that DBS of the subthalamic nucleus promotes GPi activity (Hashimoto et al. 2003). However, the level of an approximately constant inhibitory signal has relatively little impact on TC cell responsiveness to excitatory inputs, after an initial transient consisting of a few such inputs. This invariance rises because the inactivation that occurs during each TC spike and the deinactivation that occurs between TC spikes tend to balance out over the course of the transient such that the deinactivation compensates for the inactivation and allows for reliable TC responses, as long as the excitatory input frequency is not too high (Rubin and Terman 2004).

DBS leads to dispersion in TC cell failure times

The functional relevance of relay failures in TC cells will depend on how these failures are distributed across the TC population. In particular, if one TC cell bursts or fails to respond to an input but other TC cells in the population respond to this input appropriately, then the single aberrant response is unlikely to have a significant downstream effect. On the other hand, if multiple TC cells respond inappropriately to the same input, then this would be more likely to impact downstream activity.

To test the degree of temporal coincidence of TC response errors, for each case, we used a representative 3-s segment of experimental GPi recording to generate an inhibitory signal that was input to each member of a population of 40 model TC cells. All TC cells also received an identical excitatory input train, consisting of 59 pulses delivered at a frequency of 20 Hz. To enhance the realism of this computational experiment, we introduced significant heterogeneity into the TC population as described in methods. The parkinsonian and sub-therapeutic cases in these simulations are characterized by many trials in which very few TC cells achieve successful relay (Fig. 5). In contrast, in the normal and therapeutic DBS cases, there are almost no such trials (Fig. 5). More generally, the frequency distributions for numbers of TC cells achieving successful relay vary quite noticeably across the different regimes with a substantial shift in weight from trials in which most TC cells exhibit successful relay to trials in which few TC cells relay effectively and back again as GPi recording conditions switch from normal to parkinsonian without therapeutic DBS to parkinsonian with therapeutic DBS. In particular, there were statistically significant differences in the frequencies with which different numbers of TC cells responded successfully between the therapeutic DBS scenario and the other PD recording conditions (Kolmogorov-Smirnov test, P < 0.0001 for therapeutic DBS/PD as well as for normal/PD, P < 0.01 for therapeutic DBS/sub-therapeutic DBS) with a statistically insignificant difference between response frequencies in the normal and therapeutic DBS cases (P = 0.65).

FIG. 5.

TC cell failures coincide without DBS and are dispersed with DBS. A1, B1, C1, and D1: numbers of TC cells, from a heterogeneous population of 40 cells, responding successfully to each excitatory input in a train of 59 inputs (numbered 2–60, with input 1 discarded due to spurious transients), delivered at 20 Hz. For consistency, the same periodic excitatory input train was used in all cases (although we checked to ensure that qualitatively similar results held for Poisson inputs), while the GPi data used to generate the inhibition was taken either from a nonparkinsonian recording (A1), a parkinsonian recording without DBS (B1), a parkinsonian recording with sub-therapeutic DBS (C1), or a parkinsonian recording with therapeutic DBS (D1). In all cases, a successful response was defined as a response without a miss or an extra spike, as discussed in methods. A2, B2, C2, and D2: for each scenario, TC responses are collected in a histogram. To form each histogram, excitatory inputs were binned by the number of TC cells responding successfully to them. Each histogram thus shows the number of trials in which various numbers of TC cells responded successfully.

Figure 5, A2D2, summarizes this data in four histograms, one for each case. In each histogram, results are binned according to the frequency with which different numbers of TC cells responded to excitatory inputs. For example, of the 59 excitatory inputs, there were 25 inputs to which zero to eight TC cells responded successfully in the parkinsonian case without DBS (Fig. 5B2). Inspection of these plots reinforces the observation that there are many more instances of coincident TC response failures, across a large subset of the TC cell population, in the parkinsonian case in the absence of DBS than with either form of DBS, while the response failures in the presence of DBS tend to be more temporally dispersed. Moreover, this trend is a gradual one, with sub-therapeutic DBS representing an intermediate case between PD and therapeutic DBS, while the temporal dispersion of response failures in the case of therapeutic DBS resembles that of the normal case. Similar results were obtained when noise was introduced into the TC model, in addition to heterogeneity (results not shown).

Burstiness and correlation of GPi inputs both affect TC cell relay fidelity


Experimental results have shown an increase in bursting activity, as well as an increase in correlations across GPi neurons, in parkinsonian conditions, relative to normal states (Bergman et al. 1994; Brown et al. 2001; Magnin et al. 2000; Nini et al. 1995; Raz et al. 2000). However, the effect of such changes on TC relay capabilities has not been established. Because our experimental GPi data consisted of single-cell recordings, we could not use this data to assess the effect of increased correlations between GPi neurons directly. In our simulations up to this point, however, we had set gsyn for IGi→Th sufficiently high so that a single GPi input train could significantly impact TC firing. Based on this, we reasoned that the single GPi input train could be thought of as a collection of more than one, perfectly synchronized GPi signals. Correspondingly, we generated two copies of each GPi input train and divided the amplitude of the corresponding signal for each copy in half, and then we proceeded to introduce independent, normally distributed jitter into the input timing in each copy, as described in methods . We then subjected the TC cell to the jittered pair of inhibitory signals and considered how TC relay of periodic excitatory inputs varied with the amplitude of this jitter. We repeated this experiment in parkinsonian and DBS conditions, averaging over 40 jittered signals generated from a single baseline 5-s GPi data set for each case (Fig. 6, A and B).

FIG. 6.

Introducing jitter across multiple GPi signals reduces but does not eliminate the distinction between parkinsonian and DBS relay performance. Note that DBS here refers to therapeutic deep brain stimulation. A and B: the top 4 panels show GPi input signals (top traces), TC cell voltage time courses (middle traces) and excitatory inputs (bottom traces). The top 2 panels (A, 1 and 2) correspond to the DBS case, with 0 jitter on the left and σ = 0.05 on the right, while the bottom 2 (B, 1 and 2) correspond to the parkinsonian case, with 0 jitter on the left and σ = 0.05 on the right. C: error index as a function of the level of jitter amplitude σ for DBS (•) and parkinsonian (★) simulations, averaged over 40 instantiations of jitter applied to a single GPi data set for each case. D: EST vs. jitter amplitude for DBS (•) and PD (★). While EST drops with increasing jitter for both DBS and parkinsonian cases, the EST values for DBS stay well above baseline parkinsonian levels and remain at a level corresponding to significant periods of high-frequency firing, while the EST values for the parkinsonian case remain bounded away from zero. E: correlation vs. jitter amplitude for DBS (•) and parkinsonian (★) cases. Note that in the parkinsonian case, the fraction of time spent with the GPi cells simultaneously exhibiting HFE drops almost completely to 0 as jitter is increased.

The introduction of jitter within the therapeutic DBS input train had little effect on the already good TC response fidelity (Fig. 6C), although a slight smoothing of the GPi input signal, and corresponding relay enhancement, did result. Jitter amplitude did have some effect on the proportion of time during which HFE occurred in the GPi signals, as measured by the EST, and on the correlation of the pair of GPi signals in the therapeutic DBS case. However, the EST in the presence of jitter remained high (∼0.35, relative to 0.23 in the PD case without jitter), indicating that GPi inputs remained in a regime with high rates of high-frequency spiking (Fig. 6, D).

In contrast, the inclusion of jitter resulted in smoothing of the GPi input signal and, as jitter amplitude was increased, eventually yielded significant improvement in TC response fidelity in the absence of DBS (Fig. 6C). It is important to note that the introduction and gradual increase in amplitude of jitter decreased the correlation between the GPi inputs to levels near zero, but it only diminished the EST in these signals by about one third, as shown in Fig. 6, D and E, such that significant HFE remained. Indeed, the EST values for the GPi signals in the absence of DBS corresponded to bursty inhibitory time courses, featuring significant epochs with and without high-frequency spiking, for all levels of jitter. Therefore the fact that the error index dropped with increased jitter in the PD case, as can be seen in Fig. 6C, shows that input correlations likely play a role in the compromise of TC cell relay in the absence of DBS. At the same time, comparison of the PD and therapeutic DBS cases (Fig. 6, CE) shows that the error index for the PD case remains substantially above that for the DBS case, even as jitter becomes relatively large. This comparison demonstrates that the phasic or bursty nature of GPi inputs in PD, indicated here by moderate EST (Fig. 6D), also contributes significantly to the loss of TC cell relay fidelity. In summary, based on these findings, we predict that both significant correlations in GPi activity and phasic burstiness in GPi activity contribute to the compromise of TC relay fidelity in parkinsonian conditions.


To further explore the effect on TC responses induced by changes in the rate at which HFE occur and in input correlation, corresponding to the proportion of time featuring simultaneous high-frequency spiking of GPi cells, we performed simulations with computationally generated GPi input trains, for which we could control these input characteristics directly, as described in methods (Fig. 7 ). In brief, we generated two computational GPi signals, each of which depended on five stochastic spike trains, and in each spike train, HFE occurred with a rate rb. We refer to the number of spike trains that were common to both GPi signals as the number of overlaps in the simulation. For a fixed number of overlaps, we could achieve a range of correlation times by varying rb. However, with fewer overlaps, a larger rb would be required to achieve a fixed correlation level. Hence allowing different numbers of overlaps allowed us to consider more than just a one-dimensional curve in the two-dimensional space corresponding to the EST of, and the correlation between, two inhibitory input signals.

FIG. 7.

Schematic representation of the numerical generation of GPi spike times. A: each GPi cell receives and filters a combination of 5 independent random point processes, wij. An individual point process may belong to the input set of >1 GPi cell; in this example, there are 2 such overlaps, or shared wij, with w41 = w42, and w51 = w52. Varying the number of overlaps allows for control of the correlation across the inhibitory GPi inputs to the TC cell, s1 and s2 (see Eq. 3), which may affect the TC cell's responses to incoming excitatory signals. B: each wij is generated by a set of 4 Poisson processes that determine the signal's spike times and degree of burstiness or EST. Specifically, 1 process (p1) selects the times of isolated spikes, a 2nd (p2) selects the burst onset times or equivalently the times between successive bursts, a 3rd (p3) selects spike times within bursts, and a 4th (p4) selects burst durations.

We performed 3-s simulations with a range of rb values and different numbers of inhibitory input overlaps. For each simulation, we counted the number of TC misses and bad spikes (i.e., bursts or spikes not aligned with excitatory inputs, see methods) and used the results to compute the error index, according to Eq. 5, resulting from application of a 20 Hz excitatory input train. The range of error index values produced in our purely computational simulations was similar to that obtained in our simulations incorporating experimental data (compare Figs. 3 and 8), which supports the idea that our computationally generated GPi signals represent a reasonable generalization of those generated from experimental recordings. For each fixed number of overlaps, the relation between the error index and the correlation between the inhibitory inputs (achieved by varying rb) seen in our simulations is nonmonotonic: starting from small inhibitory input correlations, increases in correlations are associated with more relay errors, while starting from large correlations, further increases reduce errors in relay (Fig. 8A). A very similar trend also arises if error index is plotted against the EST of the inhibitory signals (see following text). Note also that for a fixed moderate or large value of correlation, the error index decreases as the number of overlaps decreases. For a given correlation level to occur with fewer overlaps, HFE must be present in a higher proportion of the overall inhibitory input signal; that is, the EST must be higher. Thus the cases with fewer overlaps are closer to the case of high-frequency tonic inhibition that was observed to improve relay fidelity in our other simulations (e.g., Figs. 3 and 6, therapeutic DBS case).

FIG. 8.

The error index rises and then falls again with increasing inhibitory input correlation and EST. A: error index vs. correlation, demonstrating the dependence of error index and correlation on the number of overlapping signals wij (coded by symbols and color) in the GPi input sets. Results in this and all other panels are based on simulation epochs of 3 s, with 20-Hz periodic excitation applied; similar results were obtained with Poisson input trains. B: error index vs. EST and correlation. Different symbols correspond to different numbers of overlaps (circles: 0 overlaps; triangles: 2 overlaps; squares: 4 overlaps; diamonds: 5 overlaps). The error index values are color coded such that warm colors, which occur here for moderate EST/correlation levels, correspond to high error rates and cool colors, visible here for low and high EST/correlation levels, correspond to low error rates. Note that for moderate EST, GPi firing is bursty, whereas for high EST, it is high-frequency and more tonic. C and D: the error index is decomposed into missed spikes (C) and bad spikes (D), and the dependence of each is plotted against EST and correlation. In these plots, the color bars represent the total number of occurrences observed within each 3-s simulation.

The nonmonotonic dependence of TC relay performance, measured by the error index, on correlation and EST can also be illustrated by plotting error index against both correlation and EST simultaneously (Fig. 8B). Doing so confirms that error index values peak for moderate inhibitory input EST. As noted in the preceding text, as the EST increases beyond moderate levels, the proportion of time during which high-frequency spikes are present in the inhibitory input trains increases, such that input trains approach a high-frequency, tonic spiking state (see Fig. 2) and input currents become relatively constant. In this regime, error index values decrease significantly, particularly when there are no overlaps (blue circles for large EST), which is consistent with Fig. 8A. Further, higher error rates are seen when correlations are higher, at each fixed EST, consistent with the hypothesis that synchronization of bursts of inhibition enhances their capacity to compromise relay fidelity.

Finally, we decomposed the error index into the fraction of excitatory inputs for which the TC cell fails to respond (missed spikes; see Fig. 8C) and the fraction of excitatory inputs to which the TC cell does respond but does so excessively (bad spikes; see Fig. 8D). The number of missed spikes rises significantly from low to moderate inhibitory input EST and then drops again at high EST. This number depends much more weakly on correlation, for fixed EST, than on EST itself.

Unlike missed spikes, the number of bad spikes depends strongly both on correlation and on inhibitory input EST with the highest bad spike rate occurring for relatively high correlation and moderate EST (corresponding to high burstiness). For each fixed EST, increased input correlations yield a noticeably higher rate of bad spikes. This trend makes sense because bad spikes tend to arise via a rebound effect upon the relatively abrupt withdrawal of inhibition (Rubin and Terman 2004). Such an abrupt withdrawal is more likely to occur with higher input correlations (also see Fig. 4), whereas lower input correlations lead to more smeared out input arrival times and correspondingly less abrupt changes in inhibitory currents. Similarly, for fixed correlation level, higher EST yield much lower bad spike rates, likely corresponding to the fact that with higher EST, the TC cell is subject to significant inhibition from at least one of its GPi inputs more of the time, making rebound less likely.

Taken together, the results from our computational model (Fig. 8) all support three main ideas. First, TC cell relay fidelity is compromised by inhibitory inputs that display alternations between the presence and absence of HFE with a significant correlation, or alignment of HFE, across inputs. This effect occurs through a combination of increased missed spikes and increased bad, or excessive, responses. Second, the presence of a rather tonic, high-frequency inhibitory input train, corresponding to high EST and correlation in our measures, leads to a relatively constant inhibitory current that reduces both missed and bad responses and thereby restores TC cell relay fidelity. Third, both the prevalence of the HFE and the level of the correlations in the inhibitory input structure contribute to this effect yet the contributions that these features make are distinct.


The fundamental goal of this study was to quantify how different patterns of GPi inhibition, generated from experimental recordings of normal and parkinsonian monkeys with and without DBS, affect TC relay fidelity. To this end, we subjected a Hodgkin-Huxley-type model TC cell to stereotyped excitatory signals and evaluated its ability to relay that excitatory input while under the influence of experimentally derived inhibitory pallidal modulation. We also explored a broader parameter space with computationally generated inhibitory trains in which the prevalence of high-frequency spiking episodes and the correlation structure were varied systematically. Our results show that GPi firing patterns produced in parkinsonian conditions without DBS or with sub-therapeutic DBS and, more generally, rhythmic or bursty inhibitory signals with correlations in burst timing across cells, tend to compromise the fidelity of TC cell responses to excitatory signals, relative to GPi firing patterns arising in normal conditions or in parkinsonian conditions with therapeutic DBS. More generally, improvement in TC relay fidelity was achieved by either smearing out the arrival times of correlated, bursty inhibitory signals or by converting inhibitory inputs from bursty to tonic and high-frequency. Moreover, across a model TC cell population, response failures tended to coincide temporally in parkinsonian conditions despite heterogeneity in the intrinsic characteristics of cells in the population, whereas under DBS, these failures, when they occurred, were temporally dispersed.

Multiple forms of experimental observations suggest that at least a subset of the excitatory inputs to the pallidal receiving areas of the thalamus arise from cortical areas (Guillery and Sherman 2002ac; Haber 2003). Inputs to thalamic relay cells have been classified as drivers and modulators, the former of which act on ionotropic receptors and directly induce firing and the latter of which are detectable primarily through their indirect influence on TC responses to driving signals, which may arise through action on metabotropic receptors (Sherman and Guillery 1998). Evidence has been amassed that, at least in certain thalamic areas, the excitatory drivers of thalamic relay cells represent copies of motor control signals sent from the cortex. This has led to the idea that a primary function of thalamocortical relay in the motor thalamus is to help coordinate cortical motor processing by sharing information on both motor instructions and sensory observations (Guillery and Sherman 2002b). Inhibitory inputs, on the other hand, have been posited to act as modulatory signals to TC cells (Smith and Sherman 2002). Our hypothesis about the mechanism through which parkinsonian conditions and DBS impact motor performance is consistent with this viewpoint. Specifically, our computational analysis demonstrates that differing inhibitory basal ganglia output patterns, as arise in differing nonparkinsonian and parkinsonian conditions, lead to significant differences in the ability of TC cells to relay information transmitted to these cells from other brain regions. Interestingly, we have found similar TC relay in nonparkinsonian conditions as in the parkinsonian case with therapeutic DBS. This finding suggests a way in which high-frequency stimulation of STN could restore some measure of “normal” function to the basal ganglia-thalamocortical circuit despite its profound impact on GPi activity patterns. In fact, our error index scores for the therapeutic DBS case are even lower than those based on nonparkinsonian data. We are not suggesting, however, that TC relay in isolation could be a direct measure of expected motor performance but rather that the impact of GPi firing on TC relay offers one of what are likely many mechanisms through which the effects of DBS occur. Moreover, functions that have been hypothesized to be performed by temporally precise GPi firing in normal conditions, such as termination of motor behaviors (Mink 1996; Mink and Thach 1993; Nambu et al. 2002), would presumably be disrupted by DBS, and the impacts of this disruption, as well as relevant compensation mechanisms, remain to be characterized.

The modulatory impact of GPi inhibition on TC relay in our model is mediated by the inactivation/deinactivation of spike-promoting currents, namely a sodium current (INa) and a low-threshold calcium current (IT) (see also Rubin and Terman 2004). In normal conditions, with a relatively constant inhibition from the basal ganglia to TC cells, the TC cells act in tonic mode to relay excitatory inputs, with little IT participation (Rubin and Terman 2004). In parkinsonian conditions, however, bursty inhibition leads to two effects that compromise relay, both of which are evident in Fig. 4. First, relatively abrupt rises in inhibition lead to failed relay, when excitatory inputs arrive before INa and IT can deinactivate sufficiently to overcome the inhibition. Second, subsequent deinactivation of INa and IT followed by relatively abrupt release from inhibition leads to activity bursts that do not represent excitatory input content. Finally, in therapeutic DBS conditions, although the level of inhibition to TC cells is generally higher than normal, the lack of inhibitory rhythmicity leaves IT relatively constant and therefore eliminates most rebound bursts. Moreover, the added inhibition maintains INa and IT at partially deinactivated levels, such that the added availability of these currents helps counter the direct tendency of synaptic inhibition to shunt spikes, which could otherwise lead to relay failure (Rubin and Terman 2004).

Within the literature, substantial evidence has been presented that DBS of the subthalamic nucleus (STN) suppresses or reduces somatic activity (Beurrier et al. 2001; Filali et al. 2004; Magarinos-Ascone et al. 2002; Meissner et al. 2005; Tai et al. 2003; Welter et al. 2004). Often this has been interpreted to mean that the efficacy of DBS stems from such suppression, through a removal of excessive inhibition from the targets of basal ganglia outputs (Benabid et al. 2001; Benazzouz et al. 2000; Obeso et al. 2000; Olanow and Brin 2001; Olanow et al. 2000). While this hypothesis is consistent with classical, firing-rate-based representations of information flow through the basal ganglia (Albin et al. 1989; Wichmann and DeLong 1996), it is at odds with a variety of studies showing that DBS activates areas downstream from its target site (Anderson et al. 2003; Hashimoto et al. 2003; Hershey et al. 2003; Jech et al. 2001; McIntyre et al. 2004; Miocinovic et al. 2006; Paul et al. 2000; Windels et al. 2000, 2003). From a rate-based perspective, the idea that both parkinsonian and DBS conditions lead to increased thalamic inhibition represents a paradox. This paradox may be resolved, however, by considering that DBS changes the pattern, along with the firing rate, of inhibitory inputs to thalamus (Foffani and Priori 2006; Foffani et al. 2003; Garcia et al. 2005; Meissner et al. 2005; Montgomery and Baker 2000; Rubin and Terman 2004; Terman et al. 2002; Vitek 2002). There have been some previous computational efforts to explore the details of how these varying firing patterns emerge and depend on a variety of neuronal and stimulus-related parameters (Grill et al. 2004; McIntyre et al. 2004). Building on one previous study (Rubin and Terman 2004), the work presented in this paper fills in important details of how specific changes in activity patterns induced downstream from the STN-DBS site can lead to changes in information processing through the basal ganglia-thalamocortical loop (Leblois et al. 2006; Rubchinsky et al. 2003), which would likely impact motor behavior. Interestingly, local field potential recordings from the STN of Parkinson's disease patients have shown that movement-related 300-Hz oscillations are restored by levodopa administration and contribute to related motor improvement (Foffani et al. 2003). These findings have led to the idea that high-frequency STN DBS could produce clinical benefits not only by disrupting pathological oscillations but also by driving this rhythm, at twice stimulation frequency, and thereby supporting motor processing (Foffani and Priori 2006). Our results tie in nicely with these ideas, offering one suggestion of how high-frequency oscillations in STN output could be conducive to normal information flow downstream in the network from the site of stimulation.

The incorporation of experimentally recorded GPi firing patterns into our model represents a significant advance in the computational exploration of the mechanism underlying the efficacy of DBS. As this work now stands, it represents a demonstration that in at least some subset of cells, the GPi firing pattern under parkinsonian conditions could significantly compromise TC cell relay fidelity, whereas the change in GPi firing pattern induced by therapeutic DBS could restore relay fidelity. While alterations in activity undoubtedly vary across different cells even within the same setting, the existence of changes of this type in even a subset of cells could be sufficient to affect downstream processing. One important limitation of our study, however, was the lack of simultaneous multi-unit recordings from GPi. While we were able to use computational techniques to generate simulated multi-unit inputs (Fig. 6) and to consider the impact of the experimental data on a multi-cell target population (Fig. 5), future work involving simultaneously recorded data will be performed to allow for a more direct and in-depth consideration of the activity patterns across the GPi network and the thalamic responses that these patterns induce.

Another limitation of our study was the use of a relatively simple TC cell model. We felt that it was appropriate to perform our analysis on a model that, while based on experimental data (Rubin and Terman 2004; Sohal and Huguenard 2002; Sohal et al. 2000), did not introduce undue complexity. The promising results of this study lay the groundwork for future efforts to further evaluate TC relay fidelity during DBS in more detailed, multi-compartmental TC cell models (Destexhe et al. 1998; Emri et al. 2000). In addition, the concepts considered in this work should be further explored in network models that account for the interactions of TC cells and the GPi with other brain areas, such as the thalamic reticular network (Destexhe et al. 1996; Golomb et al. 1994). In this vein, we have not assigned a specific source to the excitatory signals to the TC cell in our model, and we therefore have not considered any source-specific patterns that might be present in these signals but rather have used periodic and stochastic excitatory inputs consistent with past work (Rubin and Terman 2004). A periodic excitatory signal is nonbiological but provides the cleanest test of the effect of GPi inhibitory patterns on relay. For our stochastic excitatory inputs, we felt that, in the absence of source-specific information, it was reasonable to select the most widely used and generic form of stochastic neuronal spike train, namely a Poisson spike train. Our results do not strongly depend on input period or Poisson input rate as long as the input frequency is sufficiently high that there is little intrinsic (i.e., not input-driven) TC firing. Our results do weaken if individual excitatory input durations are made to be shorter than a few milliseconds. Thus in a situation where there is an extremely tight synchronization of excitatory inputs to a TC cell, relay fidelity differences between scenarios might be suppressed. We use a simple error index to quantify TC cell relay fidelity and a straightforward calculation of elevated spiking time to quantify the burstiness of the inhibitory signals from GPi cells. While it is possible that a more sophisticated measure would completely separate all model outputs in a single dimension, it is rather remarkable that the two simple measures used here distinguish the normal, parkinsonian, sub-therapeutic DBS, and therapeutic DBS cases so well (Fig. 3).

As noted in the preceding text, it is highly likely that cortical areas participate in driving the TC cells targeted by basal ganglia outputs. It is possible that there could be some relationship between the cortical input to these TC cells and the cortical input that enters the basal ganglia through the striatum or the subthalamic nucleus, which could then be reflected in an interdependence of the inhibitory and excitatory signals that the relevant TC cells receive. On the other hand, such a relationship might be diluted by the multi-synaptic nature of the cortico-basal-ganglia-thalamic pathway, and if DBS were applied, its effect would likely be diminished by the strong DBS signal to the STN. Given the design of our study, we did not have access to the excitatory inputs to TC cells that were present when the GPi signals were recorded experimentally. Thus in our simulations with experimental GPi signals, we could not consider the effects of any correlations between cortical signals to TC cells and cortical signals propagating through the basal ganglia, nor did we find sufficiently precise experimental characterization of such correlations to justify including them in our purely computational experiments. While it is outside the scope of this work, it would be interesting for future efforts to explore how motor signals emerging from the basal ganglia and inputs from other sources interact to shape thalamic firing patterns and how these interactions are modulated by parkinsonian rhythms. Finally, we have not analyzed the downstream responses to changes in relay across a TC cell population nor their relevance for motor performance. Future experimental work to flesh out the details of the functionally relevant inputs to, and output targets of, the thalamic cells receiving inhibition from the GPi would prove useful in pursuing such extensions of this work.


The research was partially supported by the National Science Foundation under agreement 0112050 and grants DMS0414023 and DMS0716936 to J. E. Rubin and DMS0514356 to D. Terman and by the National Institute of Neurological Disorders and Stroke Grants NS-37019 to J. L. Vitek and NS-047388 to C. C. McIntyre.


The authors thank T. Hashimoto and J. Zhang for acquiring the experimental data used in this study, G. Russo and W. Xu for processing the experimental data, and P. Hahn for assistance in selecting representative GPi recordings.


  • * Y. Guo and J. E. Rubin contributed equally to this work.

  • 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.


View Abstract