Computational Model of Recurrent Subthalamo-Pallidal Circuit for Generation of Parkinsonian Oscillations

Parkinson's disease is a movement disorder caused by dopamine depletion in the basal ganglia. Abnormally synchronized neuronal oscillations between 8 and 15 Hz in the basal ganglia are implicated in motor symptoms of Parkinson's disease. However, how these abnormal oscillations are generated and maintained in the dopamine-depleted state is unknown. Based on neural recordings in a primate model of Parkinson's disease and other experimental and computational evidence, we hypothesized that the recurrent circuit between the subthalamic nucleus (STN) and the external segment of the globus pallidus (GPe) generates and maintains parkinsonian oscillations, and that the cortical excitatory input to the STN amplifies them. To investigate this hypothesis through computer simulations, we developed a spiking neuron model of the STN-GPe circuit by incorporating electrophysiological properties of neurons and synapses. A systematic parameter search by computer simulation identified regions in the space of the intrinsic excitability of GPe neurons and synaptic strength from the GPe to the STN that reproduce normal and parkinsonian states. In the parkinsonian state, reduced firing of GPe neurons and increased GPe-STN inhibition trigger burst activities of STN neurons with strong post-inhibitory rebound excitation, which is usually subject to short-term depression. STN neuronal bursts are shaped into the 8–15 Hz, synchronous oscillations via recurrent interactions of STN and GPe neurons. Furthermore, we show that cortical excitatory input to the STN can amplify or suppress pathological STN oscillations depending on their phase and strength, predicting conditions of cortical inputs to the STN for suppressing oscillations.

Parkinson's disease is a movement disorder caused by dopamine depletion in the basal ganglia. Abnormally synchronized neuronal oscillations between 8 and 15 Hz in the basal ganglia are implicated in motor symptoms of Parkinson's disease. However, how these abnormal oscillations are generated and maintained in the dopamine-depleted state is unknown. Based on neural recordings in a primate model of Parkinson's disease and other experimental and computational evidence, we hypothesized that the recurrent circuit between the subthalamic nucleus (STN) and the external segment of the globus pallidus (GPe) generates and maintains parkinsonian oscillations, and that the cortical excitatory input to the STN amplifies them. To investigate this hypothesis through computer simulations, we developed a spiking neuron model of the STN-GPe circuit by incorporating electrophysiological properties of neurons and synapses. A systematic parameter search by computer simulation identified regions in the space of the intrinsic excitability of GPe neurons and synaptic strength from the GPe to the STN that reproduce normal and parkinsonian states. In the parkinsonian state, reduced firing of GPe neurons and increased GPe-STN inhibition trigger burst activities of STN neurons with strong post-inhibitory rebound excitation, which is usually subject to short-term depression. STN neuronal bursts are shaped into the 8-15 Hz, synchronous oscillations via recurrent interactions of STN and GPe neurons. Furthermore, we show that cortical excitatory input to the STN can amplify or suppress pathological STN oscillations depending on their phase and strength, predicting conditions of cortical inputs to the STN for suppressing oscillations.

INTRODUCTION
Parkinson's disease (PD) is characterized by progressive movement disorders caused by dopamine depletion in the basal ganglia (BG) resulting from degeneration of midbrain dopaminergic neurons. Abnormally synchronized oscillations in local field potentials (LFPs) between 8 and 30 Hz, which are called "the BG beta frequency band" (Hammond et al., 2007), are observed in the subthalamic nucleus (STN) and the internal/external segments of the globus pallidus (GPi/GPe) of PD patients (Miller and DeLong, 1987;Levy et al., 2001;Kühn et al., 2005). These abnormal oscillations are implicated in PD symptoms, such as akinesia, bradykinesia, and rigidity (Boraud et al., 2002;Brown, 2003;Gatev et al., 2006;Rivlin-Etzion et al., 2006;Hammond et al., 2007; but see Leblois et al., 2006b). We constructed a spiking neuron model of the recurrent STN-GPe circuit and investigated mechanisms for generation of parkinsonian beta band oscillations by exploring physiological parameters.
BG beta oscillations in LFPs are divided into upper-beta (>20 Hz) and lower-beta (<20 Hz) oscillations, and only the lower-beta oscillations are sensitive to dopaminergic therapy (Priori et al., 2004). Single-unit recordings from the STN in human PD patients show abnormal neuronal burst oscillations in the 10-25 Hz frequency band, which could be the source of pathological LFP oscillations (Levy et al., 2001(Levy et al., , 2002. Neuronal recordings in nonhuman primate PD models show exaggerated oscillatory burst activities in STN, GPe, and GPi that are concentrated in narrow frequency bands of 4-8 Hz and 8-15 Hz (Raz et al., 2000;Heimer et al., 2006;Wichmann and Soares, 2006;Tachibana et al., 2011), the latter of which overlaps with the frequency band recorded from human PD patients. These 8-15 Hz neuronal burst oscillations are reduced by both dopaminergic therapy and deep brain stimulation, which can alleviate PD motor symptoms Rosin et al., 2011;Tachibana et al., 2011). Recently, neuronal recordings from primate PD models under a variety of pharmacological manipulations demonstrated that STN 8-15 Hz neuronal oscillations are the origin of GPe/GPi oscillations and that they are dependent on GABAergic inhibition from GPe and glutamatergic inputs from the cortex/thalamus (Tachibana et al., 2011).
Computational modeling studies have tackled mechanisms of pathological BG oscillations (Gillies et al., 2002;Terman et al., 2002;Rubin and Terman, 2004;Leblois et al., 2006a;Hahn and McIntyre, 2010;Holgado et al., 2010;Kumar et al., 2011). However, their simulated results did not agree well with experimentally-observed neuronal properties of the BG 8-15 Hz burst oscillations (Wichmann and Soares, 2006;Tachibana et al., 2011). Furthermore, the proposed role of strengthened striatal inhibition for the induction of the BG oscillations Kumar et al., 2011) is not consistent with the experimental evidence that the blockade of GABAergic synaptic transmissions of GPe neurons from the striatum and other GPe neurons did not affect the 8-15 Hz oscillations (Tachibana et al., 2011). Therefore, the manner in which the pathological 8-15 Hz burst oscillations are generated remains unclear. Based on experimentally-observed neuronal data (Tachibana et al., 2011), we investigated mechanisms for generation and maintenance of pathological 8-15 Hz burst oscillations in recurrent STN-GPe circuits in a dopamine-depleted state. Using computer simulations of a spiking neuron model of the recurrent STN-GPe circuit, we explored the required neuronal, synaptic, and circuit mechanisms for normal firing and pathological burst oscillation. We then investigated how cortical excitatory input to the STN affects BG oscillations generated in the STN-GPe circuit and showed that for two parameter regions, the model reproduces BG neuronal activity in terms of mean firing rates, oscillation frequency, and burst discharges, corresponding to normal and parkinsonian states. The same mechanism for generation of BG pathological oscillations also predicts suppressive effects of nonoscillatory cortical inputs to the BG oscillations that explains suppression of pathological BG oscillations by movements (Amirnovin et al., 2004;Brown and Williams, 2005). We further demonstrated that oscillatory inputs from the cortex to the STN can amplify or diminish pathological BG oscillations in a phase-dependent manner, and can predict necessary conditions for their suppression.

The Model
Experimental evidence suggests that reciprocal STN-GPe interconnections and glutamatergic inputs to the STN in dopamine-depleted states are both important for generation and amplification of the parkinsonian, oscillatory burst activity of STN neurons that subsequently propagates to the GPi (Tachibana et al., 2011). Post-inhibitory rebound excitation of STN neurons is a candidate driving force of the oscillatory burst activity (Bevan et al., 2002), which could account for the importance of GPe inhibition to the STN for generation of the BG parkinsonian oscillatory burst activity of STN neurons (Tachibana et al., 2011). Short-term plasticity at excitatory synapses of GPe neurons and at GABAergic synapses of STN neurons are implicated in modulating temporal patterns of GPe and STN neurons in the normal state, and as well as in the parkinsonian state (Hanson and Jaeger, 2002;Atherton et al., 2013). Dopamine depletion in the BG is the key phenomenon for the induction of parkinsonian symptoms and BG pathological oscillations Wichmann and Soares, 2006;Tachibana et al., 2011). We considered the following neural phenomena that are related to BG dopamine depletion as candidate causes of the BG pathological oscillations: (1) reduction of GPe autonomous activity (Chan et al., 2011), (2) strengthening of GPe inhibition to the STN due to release from D2-receptor-mediated presynaptic suppression of transmitter release Johnson, 2000, 2005), and by an increase in the number of synaptic connections per GPe-STN axon terminal (Fan et al., 2012), (3) increased STN-GPe excitation to compensate for the reduction of GPe autonomous activity (Chan et al., 2011), and (4) loss of active decorrelation ability of GPe neurons (Edgerton and Jaeger, 2011). Although elevated striatal inhibition of GPe neurons by dopamine depletion (Hernandez-Lopez et al., 2000;Liang et al., 2008) is often linked to pathological BG oscillations Kumar et al., 2011), we excluded it from candidate neural phenomena related to generation of BG oscillations because experimental results showed that blockade of ionotropic GABAergic transmission either did not affect or even exaggerated GPe 8-15 Hz oscillations (Tachibana et al., 2011). Based on this evidence, we focused on a model of closed STN-GPe circuits incorporating the STN post-inhibitory rebound mechanism and dynamic synapses as mechanisms for generation of parkinsonian BG oscillations in the dopamine-depleted state. First, we explored required conditions of GPe autonomous activity and of connection strengths between the STN and GPe for generation of BG oscillations in the closed STN-GPe circuits. Later, we investigated roles of cortical excitatory inputs in the STN for BG oscillations that were generated in closed STN-GPe circuits.

Model Architecture
Anatomically, the STN and GPe are reciprocally connected. STN neurons are a major source of excitatory inputs to the GPe, and GPe neurons are a major source of inhibitory inputs to the STN. These connections are very sparse and their average connectivity is low (Kita and Kitai, 1994;Bevan et al., 1998;Sato et al., 2000;Baufreton et al., 2009). Based on these findings, we developed a random, sparsely-connected architecture of the recurrent circuit model between the STN an GPe. The architecture of our network model is shown in Figure 1A. The model consisted of 64 STN neurons and 192 GPe neurons unless otherwise noted. Each STN neuron received GABAergic projections from 12 randomlyselected GPe neurons. GPe-STN connectivity is of the same order of magnitude as experimental results (Baufreton et al., 2009). Each GPe neuron received glutamatergic projections from 6 randomly selected STN neurons. We did not include GABAergic connections between GPe neurons in our model because they were controversial: Anatomical studies showed inhibitory axonal collaterals of GPe neurons (Bevan et al., 1998;Sato et al., 2000), while functional interactions between GPe neurons were very weak or nonexistent. Each glutamatergic projection onto a GPe neuron gave a single synaptic contact, and exhibited short-term plasticity (Hanson and Jaeger, 2002; Figure 1B). Each GABAergic projection from the GPe onto an STN neuron gave N synaptic contacts that were randomly drawn from a log-normal distribution of mean = 15.50 and standard deviation = 14.53 ( Figure 1C), which fit experimental data well (Baufreton et al., 2009). Each GPe neuron had intrinsic noisy conductance (GPe Int.), as will hereinafter be described in detail. A population of cortical excitatory inputs to the STN was modeled as a stochastic spike train that obeyed a Poisson process, the mean rate of which was constant at 100 Hz, unless otherwise noted. GPe-STN GABAergic synapses exhibited short-term depression (Atherton et al., 2013). Axonal conduction delays of STN-GPe, and GPe-STN were 5 ms (Fujimoto and Kita, 1993;Kita et al., 2005).

Neuron and Synapse Models
An STN neuron was simulated by using a single-compartment, conductance-based model (Otsuka et al., 2004) that contains several sets of currents: where C m is membrane capacitance and takes 1 µF/cm 2 . V is the membrane potential. I Na is a voltage-dependent sodium current. I K is a voltage-dependent, delayed rectifier potassium current. I A is an A-type potassium current. I LCa is an L-like, long-lasting calcium current. I T is a low threshold T-type calcium current. I Ca−k is a calcium-activated potassium current. a, b, c, d1, d2, h, m, n, p, q, and r are activation or inactivation variables of these ionic currents. V Na , V K , and V Ca are the equilibrium potentials of sodium, potassium, and calcium, respectively. V leak is the resting potential of the leak current. V exc and V inh are excitatory and inhibitory synaptic reversal potentials. g Na , g K , g A , g LCa , g T , and g leak are the maximal conductances of I Na , I K , I A , I LCa , I T and I leak , respectively. Parameter values of the STN model cell were identical to the original ones except for excitatory and inhibitory synaptic reversal potentials (V exc = 0.0 mV and V inh = −84.0 mV), and the surface area of a cell, A m . Although the original model is for a unit-surface area of membrane, we set A m = 0.2 µm 2 so as to match the experimentally observed relationship between firing rate and applied current intensity (Otsuka et al., 2004). Excitatory synaptic conductance mediated by AMPAtype glutamatergic receptors, g AMPA , was modeled as the alpha function: We chose a synaptic time constant and the maximal synaptic conductance of an excitatory synapse as τ AMPA = 1.0 ms, g AMPA−Ctx−STN = 1.0 nS. Although a physiological study shows that an STN neuron receives a substantial number of excitatory inputs via NMDA-type glutamatergic receptors in addition to AMPA-type-receptor-mediated inputs (Nambu et al., 2000), we did not simulate NMDA-type conductance in this study because preliminary simulation results did not show any significant differences with or without it. The GABAergic synaptic conductance, g GABA , was modeled as the beta function: We chose τ decay = 7.7 ms and τ rise = 0.380 ms for GABAergic synapses to fit experimental data (Baufreton et al., 2009). A GPe neuron was implemented as a single-compartment conductancebased leaky integrate-and-fire model: Parameter values of GPe model neurons were as follows: C m = 0.25 pF, g leak = 16.66 nS, V leak = −70.0 mV, V exc = 0.0 mV. The absolute refractory period was set to 2.0 ms. Threshold and reset potentials were −55.0 and −60.0 mV, respectively. Synaptic conductances of GPe neurons were modeled as the beta function, after Baufreton et al. (2009): A random, sparsely-connected architecture of the recurrent circuit model between the STN and GPe. (B) Properties of three different types of short-term plasticity at STN-GPe synapses. A type of short-term plasticity is randomly chosen for a STN-GPe synapse with equal probability from these three types regardless of a post-synaptic neuron. (C) A single GPe axon gives multiple synaptic contacts on a STN neuron. Transmission probability of a unitary GPe projection is dynamically modulated by short-term depression. The number of synaptic contacts that a single GPe axon established on an STN neuron obeys the log-normal distribution.
We chose τ decay = 12.4 ms and τ rise = 5.0 ms (Hanson and Jaeger, 2002).ḡ AMPA−STN−GPe andḡ AMPA−GPe−Int represent the maximal synaptic conductances of STN-to-GPe synapses and the GPe intrinsic noisy conductance. The GPe intrinsic noisy conductance (GPe Int. for short) was introduced to model autonomous activity of GPe neurons, and was implemented as a pseudo excitatory synaptic conductance that responds to a pseudo spike train obeying a Poisson process in order to reproduce both facilitating and desynchronizing effects on GPe neuronal activity. The mean rate of a pseudo spike train was set as 100 Hz, unless otherwise noted. The dopamine depletion-related reduction of GPe autonomous activity was modeled as the reduction of the maximal conductance of the GPe intrinsic noisy conductance. Although there were other ways to implement autonomous activity, such as using singlecompartment or multi-compartment conductance-based models incorporating hyperpolarization-activated, cyclic nucleotidegated (HCN) channels, we chose this simple implementation because our focus was generation of STN parkinsonian oscillations. The means by which significant variability in intrinsic electrophysiological properties of GPe neurons (Günay et al., 2008) affects generation of parkinsonian oscillations was beyond the scope of this paper.

Short-Term Plasticity Models
Short-term depression at a GPe-to-STN GABAergic synapse was simulated using a phenomenological model: where P was the transmission probability of a unitary GPe projection. At every arrival of a presynaptic spike, P was decreased by multiplication with a decrement factor as follows: where the parameter f D < 1 and the depression of P was limited by a lower bound, P bound . This implementation smoothly increased the decrement factor from f D to 1 as the lower bound for P was reached, and smaller values of the parameter, n bound , retarded the increase of the decrement factor. At each synaptic contact of a unitary GPe projection, a spike was successfully transmitted with the probability, assuming that the transmission probability of a unitary GPe projection, P, obeyed the binomial distribution of transmissions at N contact synaptic contacts made by a unitary GPe projection. Parameters of these models were manually tuned to fit experimental data (Atherton et al., 2013; Figure 2).
We used the phenomenological model by Hanson and Jaeger (2002) to simulate short-term plasticity at STN-to-GPe glutamatergic synapses. In this model, response amplitude, A, was a product of a facilitation variable, F, and a depression variable, D (F > 1 and D < 1) that decayed exponentially: At every arrival of a presynaptic spike, the facilitation variable, F, was increased through multiplication by an incremental factor as follows: where the parameter f F > 1 and the increase of F was limited by an upper bound, F bound . At every arrival of a presynaptic spike, the depression variable, D, was decreased by multiplication by the parameter, f D (f D < 1). According to Hanson and Jaeger (2002), each synapse was randomly assigned to one of three types of plasticity, facilitation-dominant (type 1,

Data Analysis and Computer Simulations
Oscillatory activity of neurons, and oscillatory correlations of pairs of neurons were estimated by the power and cross-spectral analyses of spike trains with the shuffling method . The compensated power spectral density (PSD) was calculated by dividing the spectrum of the original spike train by the mean spectrum of locally (T = 175-225 ms) shuffled (n = 50) spike trains (Tachibana et al., 2011). A confidence level for compensated PSD was calculated based on the mean and standard deviation of the spectrum in the range of 270-300 Hz . A cell was classified as "oscillatory" when its compensated PSD had values above the confidence level (p = 0.01) in at least two adjacent frequency bins of 1,000 Hz. The compensated cross-spectral density (CSD) was calculated by dividing the CSD of the original spike train by the mean CSD of globally shuffled (n = 20) spike trains. A confidence level for compensated CSD was constructed based on the range of 270-300 Hz. A pair of neurons was classified as having "correlated oscillations" if the compensated CSD contained at least two adjacent bins within a frequency range, the values of which were above the confidence level (p = 0.01). Bursts were detected by the "Poisson surprise" method (Wichmann and Soares, 2006). A Poisson surprise value of 3 for a given sequence of at least two interspike intervals (three spikes) was used to separate burst from non-burst spike trains. All neuronal network simulations were performed in NEST (http://www.nest-initiative.org), and analyses of simulation results were performed in R and Python using Numpy, Scipy, Matplotlib, rpy2, and statsmodels libraries.

RESULTS
We constructed a spiking neuron model of the STN-GPe network by incorporating known physiological and anatomical properties of their neurons and synapses (Figures 1, 2, see Methods for details). Specifically, we incorporated short-term facilitation and depression of STN-to-GPe synapses ( Figure 1B; Hanson and Jaeger, 2002) and short-term depression of GPe-to-STN synapses (Figures 1C, 2; Atherton et al., 2013). We considered changes in strengths of STN-GPe and GPe-STN synapses and intrinsic excitability of GPe neurons (GPe Int.) as parameters possibly affected by dopamine depletion, and thereby potentially responsible for the 8-15 Hz oscillatory bursts.
We systematically searched the parameter space for regions that reproduce neuronal behaviors reported in Tachibana et al. (2011), namely, the mean firing rate, the peak power of 8-15 Hz oscillation, and the number of bursts per second (burst frequency) in normal and dopamine-depleted monkeys ( Table 1). We found that there are distinct parameter regions where simulated neuronal activities of STN and GPe exhibited properties similar to those experimentally observed in normal and parkinsonian states (Figure 3). In  (Tachibana et al., 2011). (B) Data from simulations of our model at normal and parkinsonian states shown in Figure 3.  Figure 3), where STN and GPe neuronal activities showed higher burst frequency and stronger 8-15 Hz oscillations than in the normal state. In this parameter region, mean firing rates of simulated neurons in both structures (STN and GPe) were well matched to mean firing rates of neurons experimentally recorded from dopamine-depleted, parkinsonian-state monkeys (Tachibana et al., 2011, Figure 3). The peak frequency of 8-15 Hz oscillations with maximum power ranged from 14 to 15 Hz in both structures ( Figure 4B, Table 1). This range of peak frequencies is very similar to pathological oscillations observed in dopaminedepleted monkeys (Raz et al., 2000;Heimer et al., 2006;Rosin et al., 2011;Tachibana et al., 2011). No or very few neurons exhibited significant 3-8 Hz oscillations under the conditions tested.
In the normal state, synchronous oscillatory burst discharges were not observed, and no neurons exhibited 8-15 Hz oscillations ( Figure 4A). Burst discharges were also detected in the normal state, but they occurred in a rather irregular manner and less frequently, and their durations were longer than those in the parkinsonian state. In contrast, in the parkinsonian state, neurons showed synchronous oscillatory burst discharges, with different timing in different nuclei ( Figure 4B). Each neuron exhibited 14-15 Hz burst oscillation and frequent burst firings. Cross-spectral density functions of pairs of STN and GPe neurons peaked at 14-15 Hz, as in the power spectral density functions of STN and GPe neurons ( Figure 4B). Instantaneous phase differences for STN vs. GPe pairs clustered around 90 • , indicating that the maximal firing of GPe neurons lagged behind that of STN neurons by ∼17 ms (Figure 4B, right panels), which is comparable to experimental observations in rodents (Mallet et al., 2008).

Mechanisms for Generation of 8-15 Hz Oscillatory Burst Firings
To understand neuronal and synaptic mechanisms underlying 14-15 Hz oscillatory burst firings in the parkinsonian state, we focused on a local sub-circuit in the STN-GPe network, an STN neuron (STN #46), an afferent (presynaptic) GPe neuron (GPe #41), and an efferent (postsynaptic) GPe neuron (GPe #31). Under normal conditions, STN ,and GPe neurons did not exhibit clear synchronous oscillations (Figure 5A). GPe neurons showed sustained, high-frequency discharges with occasional pauses, and the STN neuron emitted isolated spikes in an irregular manner. In spite of high firing rates of presynaptic GPe neurons, their GABAergic influences on postsynaptic STN neurons were weak ( Figure 5A) due to short-term depression of GPe-STN synapses ( Figure 5A).
In the parkinsonian state ( Figure 5B), STN neurons received strong inhibition due to burst discharges of their presynaptic GPe neurons, and were rapidly hyperpolarized so that T-currents of STN neurons were de-inactivated, as shown by an increase of the T-current inactivation variable ( Figure 5B). This induced strong post-inhibitory rebound excitation of STN neurons. Consequent burst firing of STN neurons then strongly excited their postsynaptic GPe neurons to discharge in bursts. Timing of burst discharges of efferent GPe neurons overlapped with subsequent bursts of afferent GPe neurons (Figures 4B, 5B). Although an afferent GPe neuron seldom received a reciprocal excitatory connection from its postsynaptic STN neurons under the sparse random connection architecture, these synchronous bursts of GPe neurons shut off rebound burst firing of STN neurons and then hyperpolarized them again. By repeating this cycle of events, STN-GPe circuits maintained rhythmic burst discharges.
The timing of shut-off of STN burst discharges by inhibition from GPe neurons is critical for the 14-15 Hz oscillations because post-inhibitory rebound burst discharges could Frontiers in Neuroanatomy | www.frontiersin.org last longer without active shut-offs. The slow time scale of the post-inhibitory rebound mechanism alone could lead to much slower oscillations. Although almost none of the pairs of STN and GPe neurons in the network are reciprocally connected under the connection probabilities between STN and GPe neurons in the network model, network synchronization makes presynaptic and postsynaptic GPe neurons to an STN neuron discharge in bursts in a synchronous manner.
In order to assess the role of this sparsely closed circuit of STN and GPe neurons for the generation of oscillatory burst discharges, we ran simulations of network models with different proportions of reciprocally-connected pairs of STN and GPe neurons, while the connection probability was kept constant (Figure 6A). The power of 8-15 Hz oscillations increased with an increasing proportion of reciprocally connected neurons. However, oscillations remained even when reciprocal connections were absent. We also tested whether synchronized bursts are due to highly convergent connections between STN and GPe neurons in the small network model we used. In simulations of network models with reduced connection probabilities between STN and GPe neurons by increasing the total numbers of STN and GPe neurons while maintaining the numbers of target neurons for each STN and GPe neuron constant, we observed stronger oscillatory burst activity ( Figure 6B). This indicates that generation of oscillatory burst firings is not the result of the dense connectivity between STN and GPe due to the small network size.
STN-GPe network models with higher levels of GPe intrinsic excitability, including those in the normal state, exhibited much weaker oscillations in the 8-15 Hz range. Greater GPe intrinsic excitability can reduce GABAergic conductance of STN neurons due to the strong short-term depression of GPe-STN synapses ( Figure 5B). The time constant of short-term depression was long enough to keep the synaptic transmission probability at STN GABAergic synapses almost fully depressed, indicating that short-term depression can act as a switch that turns on and off the 8-15 Hz burst oscillations in response to the level of presynaptic GPe activity. We tested the sensitivity of 8-15 Hz oscillatory burst discharges to the short-term depression parameter, n bound ( Figure 6C). Smaller n bound , which causes stronger depression, reduced the 8-15 Hz oscillations with lower levels of GPe intrinsic excitability, and larger n bound enhanced the oscillations.
Another possible effect of higher GPe intrinsic excitability is enhanced stochasticity in GPe firing, because we implemented that using noisy conductance that obeyed a Poisson process (see Method). To address this issue, we investigated the effect of varying GPe intrinsic noisy conductance by increasing the mean rate of the Poisson process and reducing the maximal conductance inversely so as to keep the mean value constant. With decreased variance of GPe noisy intrinsic conductance, a slight enhancement of 8-15 Hz oscillations occurred at intermediate GPe intrinsic conductances (0.75 nS), but very little enhancement was observed at higher levels (1.5 nS), even when the conductance did not fluctuate stochastically, but remained constant at the mean value ( Figure 6D). This result suggests that stochasticity in GPe intrinsic conductance is not a major factor in preventing the occurrence of 8-15 Hz oscillations.

Roles of Cortical Excitatory Inputs on STN Parkinsonian Oscillations
Next, we investigated how cortical excitatory inputs to the STN affect parkinsonian oscillations. Anatomically, layer-V pyramidal-tract neurons in the motor cortex are a major excitatory source for the STN Kita and Kita, 2012). Layer-V pyramidal-tract neurons exhibit spontaneous firing and phasically increase their activities just before and during movements (Turner and DeLong, 2000;Pasquereau and Turner, 2011). First, we investigated the effects of cortical tonic background excitatory inputs to STN neurons on STN parkinsonian oscillations. We used homogeneous Poisson spike trains of fixed mean rate (100 Hz) as the background excitatory inputs to STN neurons, which were not correlated with BG parkinsonian oscillations. In contrast to experimental evidence (Tachibana et al., 2011), reducing synaptic strengths of background inputs to the STN enhanced 8-15 Hz neuronal oscillations in the STN ( Figure 7A) and GPe (data not shown). This is because tonic excitatory inputs to STN neurons antagonized GPe inhibition, which depolarized STN neurons and diminished post-inhibitory rebound excitation. The suppressive effect of the background excitatory inputs on the parkinsonian oscillations suggests that stronger cortical excitatory inputs can reduce BG parkinsonian oscillations. This is consistent with experimental reports showing that pathological oscillations in LFPs recorded from the STN of PD patients were reduced by voluntary movements (Amirnovin et al., 2004;Brown and Williams, 2005) because just before and during voluntary movement, cortical inputs to the STN are phasically enhanced (Turner and DeLong, 2000).
Simulation results indeed showed that parkinsonian neuronal oscillations dropped at about the time of phasic increases (up to 1000 Hz) of excitatory inputs to STN neurons (Figures 7B,C). Reduction of STN parkinsonian neuronal oscillations correlated with reduction of STN T-current conductances that underlie post-inhibitory rebound excitation (r 2 = 0.26, p < 0.001; Figure 7D). Smaller phasic increases (up to 250 Hz or 500 Hz) of excitatory inputs to STN neurons also showed drops of STN oscillations (Friedman test, p < 0.001; pairwise comparisons using Wilcoxon rank sum test, p < 0.001, data not shown).
We also investigated effects of oscillatory cortical inputs to STN neurons on STN parkinsonian oscillations. Because the manner in which cortical oscillations relate to STN parkinsonian oscillations is still unclear (Brittain and Brown, 2014), we simply assumed that STN neuronal oscillations propagate through the STN-GPi-thalamus-cortex loop and go back to the STN, while maintaining their oscillation frequency. We investigated the effect of the time delay in the STN-GPi-thalamus-cortex feedback loop on STN 8-15 Hz oscillations by introducing delay lines via which STN neurons were recurrently connected. More precisely, each STN neuron was randomly connected with 16 afferent STN neurons (average) and spikes were stochastically transmitted (p = 0.25) via delay lines. Synaptic strength of delayed feedback inputs was the same as that of background inputs to the STN. Simulation results showed that STN 8-15 Hz neuronal oscillations were affected by both the delay and the strength of the recurrent excitatory input to STN neurons (Welch's one-way ANOVA F = 4741.3, df = 17, p < 0.001). When the delay was within 60-80 ms, which was close to the cycle of STN parkinsonian oscillations, the STN oscillation in the parkinsonian state was stronger with feedback recurrent inputs than without (Figure 7E, pair-wise t-test with Bonferroni adjustment, p < 0.001, n = 512). In contrast, when the delay was shorter than 60 ms or longer than 80 ms, feedback recurrent inputs diminished STN oscillations to less than those without feedback. The power of STN oscillations in the parkinsonian state with feedback recurrent inputs, and with a delay between 60 and 80 ms was even greater than with a blockade of AMPA transmission both at background and at recurrent inputs to the STN (Figure 7E, pair-wise t-test with Bonferroni adjustment, p < 0.001, n = 512). Feedback excitatory inputs did not FIGURE 7 | Continued the STN neurons depended on different GPi-thalamo-cortical delays of cortico-STN feedback oscillatory inputs. Each statistical value was calculated from pooled spike trains of STN neurons in 8 independent simulations (total n = 512 per condition). Red hatched and blue bars correspond to network parameters for the parkinsonian and normal conditions shown in Figure 3, respectively, except as otherwise noted. "delay" indicates the GPi-thalamus-cortex delays of AMPA-mediated recurrent feedback oscillatory inputs to the STN. Presence and absence of tonic background (bkg) and recurrent (rec) cortical excitatory inputs are shown as o and -, respectively. The leftmost bar represents data from the parkinsonian state without cortical recurrent inputs to the STN, shown in Figure 3, and the second bar from the left represents data of the parkinsonian state without either the tonic background or recurrent cortical inputs to the STN. Asterisks, hash marks, and plus marks represent significant differences (pair-wise t-test with Bonferroni adjustment, p < 0.001, n = 1,024) from the parkinsonian state without AMPA-mediated, recurrent oscillatory inputs to the STN, from the parkinsonian state without either the AMPA-mediated, background or recurrent oscillatory inputs to the STN, and from the parkinsonian state with AMPA-mediated, background inputs to the STN, and recurrent oscillatory inputs to the STN with a delay of 70 ms, respectively.
affect STN 8-15 Hz oscillations in the normal state ( Figure 7E, pair-wise t-test with Bonferroni adjustment, p = 1.0, n = 512). These results indicate that cortical excitatory inputs affect STN parkinsonian oscillations in a complex manner: cortical non-oscillatory inputs are suppressive to STN parkinsonian oscillations, but feedback oscillatory inputs can amplify or diminish STN parkinsonian oscillations in a delay-dependent manner. The complex effects of cortical excitatory inputs on STN parkinsonian oscillations explain two lines of experimental evidence suggesting suppressive or facilitative effects of cortical-STN excitatory inputs, respectively (Amirnovin et al., 2004;Brown and Williams, 2005;Tachibana et al., 2011).

DISCUSSION
In this study, we explored mechanisms of generation of parkinsonian 8-15 Hz neuronal burst oscillations in the STN-GPe circuit by building a spiking neural network model, incorporating anatomical and physiological constraints. Through a systematic search of the parameter space, we identified parameter regions that reproduce normal and parkinsonian states as reported by Tachibana et al. (2011). The parkinsonian state was characterized by reduced intrinsic excitability of GPe neurons and increased conductance of GPe-STN synapses. Reduction of autonomous activity of GPe neurons by dopamine depletion paradoxically elevated inhibition of STN neurons via strong short-term depression and consequently induced strong post-inhibitory rebound burst firing of STN neurons, which were then shaped into synchronous oscillations in the 14-15 Hz frequency range via interactions between STN and GPe neurons. Characteristics of neuronal activities observed in our simulations agreed well with experimental findings. First, generation of neuronal burst oscillations is critically dependent on GABAergic inhibition of STN neurons (Tachibana et al., 2011) and is sensitive to a level of GABAergic inhibition of STN neurons (Figure 3). Second, the frequency range of these neuronal burst oscillations is in the 8-15 Hz frequency band (Figure 4B), in which pathological oscillations were observed in the BG of dopamine-depleted monkeys Wichmann and Soares, 2006;Tachibana et al., 2011). Third, mean firing rates of neurons in the model exhibited in simulations under normal or dopamine-depleted states (Figure 3) accord well with those experimentally observed in the BG of monkeys (Tachibana et al., 2011). Fourth, inactivation of the STN diminishes 8-15 Hz neuronal oscillations in the GPe and GPi (Tachibana et al., 2011), which is consistent with our simulations with smaller synaptic strengths of STN projections (Figure 6). Fifth, acute administration of the dopamine precursor, L-DOPA, reduces the power of 8-15 Hz oscillations, the strength of burst discharges in the STN and GPi, and firing rates in the STN (Tachibana et al., 2011). These experimental results are consistent with our simulation results (Figures 3, 4) if L-DOPA reduces GABAergic inhibition of the STN by restoring D2-receptormediated presynaptic suppression Johnson, 2000, 2005), but does not acutely elevate GPe autonomous activity through dopamine-dependent expression of HCN channels, which is supposed to be a slower process (Chan et al., 2011). Sixth, burst frequency in the STN and GPe (but not in the GPi) in the model were higher in the pathological state than that in the normal state (Figures 3, 4), which is consistent with experiments (Wichmann and Soares, 2006;Tachibana et al., 2011). Seventh, our simulation results suggest that contributions of cortical excitatory inputs to the STN for PD burst oscillations are complex and depend on their temporal patterns and firing rates (Figure 7). Enhancing non-oscillatory cortical inputs suppressed STN PD oscillations, which is consistent with experimental evidence from PD patients (Amirnovin et al., 2004;Brown and Williams, 2005), although there are also experimental results showing that the blockade of ionotropic glutamatergic transmission in the STN diminished STN PD oscillations (Tachibana et al., 2011). Additional oscillatory inputs, the frequency of which was similar to STN PD oscillations, could amplify or diminish STN PD oscillations in a phase-dependent manner. These amplifying effects of oscillatory cortical inputs are consistent with experimental results (Tachibana et al., 2011).
We verified the critical role of strong post-inhibitory rebound excitation of an STN neuron for generation and maintenance of PD 8-15 Hz neuronal burst oscillations (Figure 5). The same mechanism was already explored in the pioneering study of a spiking neuron model for PD oscillations , which explained the slower oscillations, but not the 8-15 Hz oscillations. The slowly-changing nature of the rebound mechanism is ascribed to the slower oscillations (Holgado et al., 2010). In our model, in response to STN rebound firing, GPe neurons actively terminate the rebound excitation of the STN despite sparsely closed loops between STN and GPe neurons, where each population was synchronously activated. This network mechanism enabled faster neuronal oscillations in spite of the slow rebound mechanism. The suppressive role of excitatory inputs on the STN for the 8-15 Hz oscillations we observed can be explained as the damping of the post-inhibitory rebound excitation by preventing required hyperpolarization. It is consistent with peri-movement, phasic activation of layer-V pyramidal tract neurons in the motor cortex projecting to both the spinal cord and the STN (Turner and DeLong, 2000). It is also consistent with the therapeutic effects of deep brain stimulation to the STN that is able to suppress 8-15 Hz oscillations by virtue of its high-frequency, direct stimulation of cortical afferents to the STN (Miocinovic et al., 2006).
The decline of GPe autonomous activity due to down regulation of the HCN channel caused by dopamine depletion (Chan et al., 2011) is the key phenomenon in our model for converting normal neuronal activity into a pathological state. Normal autonomous activity of GPe neurons could have two possible effects: maintenance of a high level of neuronal firing, and decorrelation of GPe neuronal activities by collateral inhibition (Edgerton and Jaeger, 2011). Our simulations showed that GPe neuronal firing maintained at a high level reduces GABAergic inhibition of the STN via strong short-term depression (Atherton et al., 2013) and prevents the generation of pathological oscillations. Once GPe neuronal firing is reduced due to reduced autonomous activity caused by dopamine depletion, GABAergic inhibition becomes less depressed and induces the rebound excitation of STN neurons, and consequently pathological oscillations occur. Decorrelated GPe neuron activities may prevent the generation or maintenance of 8-15 Hz burst oscillations, but it showed limited effects in a PD state in our simulations. Decorrelation of GPe neurons might be more implicated in normal BG functioning.
Although elevated striatal inhibition of GPe neurons by dopamine depletion (Liang et al., 2008) has been considered the key cause of PD pathologies (Albin et al., 1989;DeLong and Wichmann, 2007), we did not explore its role relative to pathological oscillations because of negative experimental evidence showing that blockade of ionotropic GABAergic transmission in the GPe does not change or elevate the power of the pathological oscillations (Tachibana et al., 2011). However, elevated striatal inhibition could have a similar effect, since the decline of GPe autonomous activity in our model reduced GPe firing to cause pathological oscillations.
Our simulation results showed that in-phase oscillatory excitatory inputs from the cortex amplify pathological oscillations in the STN, as observed experimentally (Tachibana et al., 2011). Indeed, cortical and STN oscillatory activities are mainly coherent in the low (12-20 Hz) and high (20-35 Hz) beta range (Hirschmann et al., 2011). Whether STN PD lower-beta oscillations proliferate throughout the cortex-BG-thalamus loops as an entrainment loop is currently a matter of debate (Brittain and Brown, 2014). However, enhancement of cortical lower-beta oscillations in PD patients (Pollok et al., 2012) and neuronal lower-beta (and upper-beta) oscillations of pyramidal tract neurons in the motor cortex of dopamine-depleted primates (Pasquereau and Turner, 2011) are consistent with the entrainment loop hypothesis. The required delay from the STN to the cortex might be too large considering just conduction and synaptic delays. However, because GPi projections to the thalamus are inhibitory, oscillations could be delayed by more than the sum of conduction and synaptic delays, as observed in closed STN-GPe circuits. The reduction of STN PD lower-beta oscillations by out-of-phase cortical oscillations in our simulations is a result of damping of STN post-inhibitory rebound excitation by counteraction of phasically-elevated excitatory inputs from the motor cortex to inhibition from the GPe. It is similar to the suppression by dynamically increased cortical inputs during movements (Turner and DeLong, 2000), and account for mechanisms of closed-loop deep brain stimulation (Rosin et al., 2011).
Our model has several limitations in reproducing experimentally observed phenomena. First, our model did not show neuronal oscillations in the 3-8 Hz frequency band in the BG of dopamine-depleted monkeys Wichmann and Soares, 2006;Tachibana et al., 2011), although the relationship of such oscillations to dopamine depletion is controversial. Second, burst frequency of the STN in the parkinsonian state in our simulations was lower than in experimental observations (Wichmann and Soares, 2006;Tachibana et al., 2011).
In summary, our model provides plausible mechanisms for generation and maintenance of PD oscillations in the BG. Our model also showed that cortical inputs to the STN can amplify PD oscillations and it predicted required conditions for amplification or reduction of pathological BG oscillations.

AUTHOR CONTRIBUTIONS
OS and KD planned the study. YT and AN provided physiological data. OS constructed the model and performed numerical simulation. All authors discussed the results and OS and KD wrote the manuscript.