Autonomous oscillations and phase-locking in a biophysically detailed model of the STN-GPe network

The aim of this study was to understand the relative role of autonomous oscillations and patterning by exogenous oscillatory inputs in the generation of pathological oscillatory activity within the subthalamic nucleus (STN) - external globus pallidus (GPe) network in Parkinson’s disease. A biophysically detailed model that accounts for the integration of synaptic currents and their interaction with intrinsic membrane currents in dendritic structures within the STN and GPe was developed. The model was used to investigate the development of beta-band synchrony and bursting within the STN-GPe network by changing the balance of excitation and inhibition in both nuclei, and by adding exogenous oscillatory inputs with varying phase relationships through the hyperdirect cortico-subthalamic and indirect striato-pallidal pathways. The model showed an intrinsic susceptibility to beta-band oscillations that was manifest in weak autonomously generated oscillations within the STN-GPe network and in selective amplification of exogenous beta-band synaptic inputs near the network’s endogenous oscillation frequency. The resonant oscillation frequency was determined by the net level of excitatory drive in the loop. Intrinsically generated oscillations were too weak to support a pacemaker role for the STN-GPe network, however, they were considerably amplified by sparse cortical beta inputs when their frequency range overlapped and were further amplified by striatal beta inputs that promoted anti-phase firing of the cortex and GPe, resulting in maximum transient inhibition of STN neurons. The model elucidates a mechanism of cortical patterning of the STN-GPe network through feedback inhibition whereby intrinsic susceptibility to beta-band oscillations can lead to phase locked spiking under parkinsonian conditions. These results point to resonance of endogenous oscillations with exogenous patterning of the STN-GPe network as a mechanism of pathological synchronization, and a role for the pallido-striatal feedback loop in amplifying beta oscillations. Author summary Exaggerated beta-frequency neuronal synchrony is observed throughout the basal ganglia in Parkinson’s disease and is reduced with medication and during deep brain stimulation. The power of beta-band oscillations is increasingly used as a biomarker to guide antiparkinsonian therapies. Despite their importance as a clinical target, the mechanisms by which pathological beta-band oscillations are generated are not yet clearly understood. In vitro electrophysiological recordings support a theory of enhanced phase locking of the reciprocally connected subthalamo-pallidal network to beta-band cortical inputs but this has not yet been clearly demonstrated in a model. We present a new model of the subthalamo-pallidal network consisting of biophysically detailed cell models that captures the interaction between synaptic and intrinsic currents in dendritic structures. The model shows how phase locking of subthalamic and pallidal neurons and exaggerated bursting in subthalamic neurons can arise from the interaction of these currents when the balance of excitation and inhibition is changed and how phase locking is amplified under specific phase relationships between cortical and striatal beta inputs.

Pathological oscillations in the basal ganglia-thalamocortical (BGTC) network have 2 long been implicated in the motor symptoms of Parkinson's disease. Beta-band (15-30 3 Hz) oscillations are strengthened consistently after dopamine depletion (DD) both in 4 Parkinson's disease (PD) patients and parkinsonian animal models [1][2][3], and are 5 reduced by deep brain stimulation and pharmacological interventions that alleviate 6 parkinsonian motor symptoms [4][5][6][7]. The magnitude of beta oscillations has also been 7 shown to correlate with the severity and degree of improvement of bradykinetic/akinetic 8 motor symptoms and rigidity [4,8]. 9 The reciprocally connected subthalamo-pallidal network is a key site in the basal 10 ganglia (BG) in which beta-band oscillations are manifest in Parkinson's disease [2,9]. 11 It is not clear whether it plays an active part in generating beta-band oscillations, nor 12 whether it amplifies or merely sustains them. Neither is it fully understood how 13 beta-band oscillations relate to other pathological patterns of neural activity in the 14 subthalamic nucleus (STN) and external globus pallidus (GPe) that correlate more 15 strongly with parkinsonian motor symptoms, notably increased neural bursting [10,11]. 16 It is clear, however, that interventions in the loop and its afferents that reduce 17 beta-band oscillations [12] or bursting [13][14][15], in the form of precise chemical lesions or 18 optogenetic stimulation targeted at the subcellular level, lead to improvements in motor 19 symptoms. Similarly, the STN and GPe are effective targets for deep brain stimulation 20 (DBS) [16,17], which likely activates a more diverse set of neuronal elements [18]. 21 Although at present it is unknown whether a causal relationship exists between 22 beta-frequency activity and parkinsonian motor symptoms, due to the strong 23 correlation between the two it is important to understand the biophysical mechanisms 24 underlying these pathological activity patterns. 25 Computational models provide a valuable tool with which to explore different 26 hypotheses regarding the origin of pathological synchronization in the BG and their 27 suppression by means of therapeutic interventions. The most prominent hypothesis 28 generally emphasize the importance of a particular feedback loop between 29 subpopulations of cells within the multiple interconnected loops that constitute the full 30 BGTC network. Different models have placed the origin of beta and sub-beta band 31 oscillations in the STN-GPe network [19][20][21][22], in striatal or pallidostriatal 32 circuits [23,24], or in the full BGTC loop [25][26][27][28] (reviewed in [29,30]). These models, 33 based on experimentally observed oscillations in different subcircuits of the BG, show 34 that many of them are prone to oscillate, be it through intrinsic pacemaking or 35 susceptibility to an extrinsic rhythm. 36 The STN-GPe network was an early focus of modelling studies due to its reciprocally 37 connected structure and because it can generate low frequency oscillations 48 Cortical patterning of the STN-GPe network by means of feedback inhibition 49 provides a proposed mechanism for this functional coupling [9,12,[39][40][41]. In this 50 hypothesis, weak oscillatory activity arriving via cortico-STN afferents is amplified in 51 the STN-GPe network when feedback inhibition by GPe is offset in phase to cortical 52 excitation. While such feedback-mediated oscillations have been observed in 53 vivo [39,42], the ability of the network to generate autonomous oscillations and its 54 resonant response properties are still poorly understood. 55 To explain oscillatory properties, previous modeling studies have focused on 56 alterations in connection patterns and strengths within or between nuclei, typically 57 represented by mean-field or single-compartment spiking neuron models. While such 58 models are computationally efficient, they may not fully capture the role of intrinsic 59 properties of neurons in shaping pathological activity patterns (reviewed in [30]). 60 Although cell-specific ion channels can be used, single-compartment neuron models 61 lump together ion channels and synapses in one isopotential compartment in a way that 62 may not capture the complex dynamics that arise when non-uniformly distributed ion 63 channels [43] interact with synapses associated with distinct subcellular 64 regions [15,44,45]. Hence they may not fully account for the mechanisms contributing 65 to cortical patterning and the role that synaptic-ionic current interactions play in 66 sustaining beta-band oscillations and excessive burst firing. 67 Recent evidence suggests that following dopamine depletion the balance of excitatory 68 and inhibitory synaptic currents in STN neurons is shifted toward inhibition [46,47] 69 which is known to promote burst responses through increasing the availability of Ca 2+ 70 and N a + channels that are de-inactivated at hyperpolarized potentials [39]. In the GPe 71 an increase in inhibition, caused mainly by strengthening of striato-pallidal afferents, is 72 also believed to play a role in generating pathological oscillations as demonstrated in 73 model simulation [19,21,32,48]. Based on experimental data, increased GPe inhibition 74 has been suggested to cause increased engagement of HCN channels [49], which are 75 involved in phase resetting and controlling the regularity of firing. However, whether 76 the functional coupling between BG nuclei is also moderated by the 77 excitation-inhibition balance is not fully understood. 78 The aim of this study was, therefore, to understand the conditions under which 79 cortical patterning of the STN-GPe network arises in a biophysically plausible manner 80 and how this is related to the balance of excitation and inhibition in each population. 81 The STN-GPe network was modeled using biophysically detailed multi-compartmental 82 cell models of STN and GPe neurons that capture the interaction between synaptic and 83 intrinsic currents distributed within the dendritic structure and involved in autonomous 84 pacemaking and bursting [43,50] (Fig. 1). The generation of oscillations both 85 autonomously within the network and in response to beta frequency inputs from the 86 cortex (CTX) and striatal medium spiny neurons (MSN) was examined as the balance 87 of excitation and inhibition within the network was systematically varied, and 88 oscillatory inputs with varying phase relationships were added. A better understanding 89 of the relative contribution of these different factors and their interaction has the 90 potential to improve understanding of the mechanism of action of existing 91 anti-parkinsonian therapies, including deep brain stimulation (DBS) and to guide the  were modeled using multi-compartmental neuron models. Cortical projection neurons (CTX) and striatal medium spiny neurons (MSN) were modeled as spike generators. B: Branching structure of the STN neuron model and representative synapses by afferent type, indicating subcellular distribution of synapses. Cortical glutamergic afferents synapse primarily onto thin dendrites, distally relative to the soma, but NMDA receptors with faster NR2A subunits mainly target the soma and proximal areas. Pallidal GABAergic afferents target proximal areas of the cell. C: Branching structure of the GPe neuron model and representative synapses by afferent type. GABAergic GPe-GPe collaterals mainly target somata and proximal dendrites, whereas glutamergic afferents were placed in distal regions. Full details of the model are provided in the Methods.

94
The STN-GPe pacemaker hypothesis was first investigated by modeling cortical inputs 95 to the STN as Poisson spike generators without any periodic or oscillatory component projections between nuclei. The ratio of total excitatory to inhibitory synaptic currents 103 (E/I ratio) was shifted by scaling the peak conductance of all synapses belonging to a 104 given projection known to be strengthened or weakened by dopamine depletion. The 105 role of additional oscillatory inputs entering the STN-GPe network via the indirect 106 striato-pallidal pathway and their phase relationship to cortical inputs was then 107 investigated.

108
The excitation-inhibition balance in the STN sets the oscillation 109 frequency of the STN-GPe network and firing mode of STN 110 neurons.

111
Sweeping the synaptic conductance of cortico-STN afferents moved the STN-GPe 112 network through a diverse range of firing patterns (Fig 2). Increasing the synaptic 113 conductance of cortical Poisson-distributed inputs caused a proportional increase in net 114 excitatory current to the STN resulting in increased excitatory input from STN to GPe. 115 However, the GPe population exhibited a saturating population firing rate curve due to 116 two negative feedback mechanisms that counteract the effect of increasing STN 117 excitation and have a homeostatic effect on the GPE's ratio of excitation to inhibition: 118 short-term depression of STN to GPe synapses [51] and reciprocal inhibition through 119 intra-GPe colaterals. As a result, the increase in feedback inhibition from the GPe did 120 not keep pace with the rise in excitation from the cortex (Fig 2.Ai-ii). At the lower end 121 of the synaptic conductance scale, both excitation and inhibition were low and the ratio 122 of excitatory to inhibitory synaptic currents (E/I ratio) was balanced towards inhibition 123 (E/I < 1). STN and GPe neurons showed high spike train synchronization at 12 Hz

129
The synchronization frequency of GPe neurons increased across the parameter sweep 130 but oscillation power was low and the oscillation was weakly transmitted to STN 131 neurons in the baseline model (Fig 2.Bi-Bii). Although synchronization of STN neurons 132 was weak throughout most of the parameter range, spikes showed a tendency to occur 133 in a fixed phase relationship with GPe as evidenced by the alignment of individual and 134 population phase vectors (Fig 2.Cii,Dii,Eii) and by the phase histogram (Fig 5.Bi).

135
Only when both the excitatory and inhibitory drive to STN were low, and there was a 136 net inhibition to STN neurons, did the STN and GPe synchronize strongly (Fig 2.Bi-Bii, 137 gmax scale = 0.2; Fig 2.Ci-Ciii). This was because in the baseline model, GPe to STN 138 inhibition was dominated by slow GABA B receptor (GABA B R) mediated synaptic 139 currents that promote low-frequency bursting and counteract faster synchronization in 140 STN. These slow GABA B R-mediated currents are facilitated by faster pre-synaptic 141 spiking in a near-sigmoidal fashion and thus tend to dominate at higher GPe firing rates. 142 When GPe to STN inhibition was shifted to the faster GABA A synapses, low-frequency 143 bursting was reduced, beta-band synchronization was strengthened and the frequency inhibition because of the negative feedback mechanisms in the GPe described above.

148
Excitation and inhibition were relatively balanced (0.7 ≤ E/I ≤ 1.1) but the total 149 inhibitory current was higher. As detailed above, higher inhibition in the dendrites 150 leads to higher availability of dendritic ion channels underlying burst responses. Hence, 151 the shift to higher inhibition brought STN neurons into a slow burst firing mode with 152 sparse, strong bursts with higher intraburst firing rate (Fig 2.Dii). These low-frequency 153 fluctuations in firing rate were transmitted to GPe neurons and showed up clearly in the 154 power spectrum of both nuclei (Fig 2.Bi-Bii). As the CTX to STN conductance was 155 increased further, resulting in a net excitation to STN, STN firing rates increased 156 accordingly and became more uniform, reflected in a lower coefficient of variation and a 157 convergence of the background and intra-burst firing rates (Fig 2.Ei). Increasing the synaptic conductance of the inhibitory GPe-GPe projection moved 159 the STN-GPe network through a similar trajectory of firing patterns while shifting the 160 excitation-inhibition balance in STN and GPe in opposing directions (Fig. 3). The 161 strong low-frequency burst firing mode of STN cells was engaged within the same 162 interval of E/I ratios (Fig. 3, Ai, Bi, F), and firing became more temporally uniform as 163 this ratio increased to favour excitation ( Fig. 3, Ci, Di, Eii, F). As with the parameter 164 sweep of the CTX-STN conductance, the mean firing rate of the GPe population was 165 less sensitive to the synaptic conductance than that of the STN population. The 166 negative feedback structure of the STN-GPe network was engaged to maintain the mean 167 GPe firing rate within a narrow range: the increase in collateral GPe-GPe inhibition 168 was offset by the increased excitation by STN cells as they were disinhibited (Fig. 3.Aii). 169 GPe neurons again synchronized more strongly than STN neurons: the 13-30 Hz peak 170 was prominent in its power spectral density (PSD) (Fig. 3.B) and single unit spikes were 171 more tightly locked to the phase of the oscillation compared to STN (this was true 172 whether the instantaneous phase was extracted from the STN or GPe population, data 173 not shown). Although STN cells were overall more likely to spike during a specific phase 174 of the oscillation, their tendency to emit longer bursts rather than short bursts or  Low-frequency bursting in STN neurons is engaged for weak GPe-GPe inhibition leading to strong GPe to STN inhibition (Bi). GPe neurons synchronized at beta frequencies (Bii). C-D: Representative spike trains and phase vectors for STN (column i, green) and GPe population (column iii, red) for two values of the GPe to GPe conductance (scale 0.33; 2.0 in rows C; D respectively). Column ii shows phase vectors of the STN and GPe populations (in green; red, respectively, mean population vectors plotted as thick solid lines and cell vectors as thin transparent lines) reflecting phase locking to the instantaneous GPe phase. STN neurons show a preferred phase of firing but weak beta synchronization, while GPe neurons synchronize more strongly. Ei: Population vector length and angle of STN and GPe population (green; red, respectively). F: Metrics that characterize bursting in STN neurons: median burst rate, intra-burst firing rate, and coefficient of variation of ISIs across all STN cells. High STN inhibition favors strong bursts (high intra-burst firing rate) and lower inhibition regularization the firing rate (lower CV). Increasing the synaptic conductance of the GPe to STN projection shifted the firing 183 mode of STN neurons towards sparse, longer bursts with high intra-burst firing rate set 184 against a lower background firing rate, characterized by a high coefficient of variation of 185 inter-spike intervals (ISIs) (Fig 4.F). Bursting was periodic at low frequencies (˜2-5 Hz) 186 but was not synchronized between cells. While STN neurons were weakly entrained to 187 the beta oscillation in the GPe they preferentially fired in an interval leading the GPe 188 by approximately 65 degrees (Fig 4.E). Longer bursts with increased firing rate were 189 promoted by higher availability of voltage-sensitive N a + and Ca 2+ channels under 190 more hyperpolarized conditions [39,43,52]. This shift towards low-frequency strong 191 bursting was accompanied by increased sub-beta band power, and intra-population  Endogenous oscillations in the STN-GPe network are strengthened by shifting the GPe to STN synaptic current from slow GABA B receptors to fast GABA A receptors. A: Reduced GABA B R-mediated currents and increased GABA A R-mediated currents lead to strengthening of endogenous oscillations. The GABA B conductance of the GPe to STN projection was decreased by 50% and the GABA A conductance was increased progressively (horizontal axis). B-D: Effect of shifting the inhibitory current from GABA B to GABA A receptors while keeping the E/I ratio in each population approximately constant. The GABA A conductance was increased by a factor six and the GABA B conductance was adjusted so that the E/I ratio was close to that in the baseline model (left column: baseline model, E/I ratio was 0.89; 0.97 in STN, GPe respectively; right column: scaled conductances, ratio was 0.89; 0.93). B: Phase histograms of STN (green) and GPe (red) neurons in baseline model (left; Bi-Bii) and model with scaled conductances (right; Biii-Biv). C-D: Representative spike trains of STN (green) and GPe (red) neurons in baseline model (left; Ci-Di) and model with scaled conductances (right; Cii,Dii). E: Repeat of the parameter sweep of the CTX to STN conductance in Fig. 2, with scaled GABA conductances. The GABA B conductance of GPe to STN synapses was halved, and the GABA A conductance was doubled. Mean PSD of somatic membrane voltages of STN (Ei) and GPe (Eii) neurons. The oscillation is stronger compared to the baseline model, is reliably transmitted to STN, and the frequency is shifted by the level of cortical excitation.

STN-GPe network shows resonant properties and phase locks to 211 cortical beta inputs 212
According to the cortical patterning hypotheses, coherent beta oscillations in the STN 213 and GPe arise due to phase locking of the STN-GPe network to cortical beta-range 214 inputs that are amplified within the loop. To assess the degree of phase locking to such 215 cortical rhythms and its dependence on intrinsic parameters of the STN-GPe network, 216 the network was simulated with cortical inputs modeled as oscillatory beta bursting  Striatal microcircuits exhibit beta-band oscillations in healthy primates [54] and 237 parkinsonian mice models [23] and have been hypothesized to be part of the pacemaking 238 April 11, 2019 8/33 . The x-axis shows the estimated balance of excitatory and inhibitory currents measured from 3 recorded neurons in each population. Maximum phase locking to cortical frequency occurs at approximately the same E/I ratio of maximum autonomous power at the same frequency.
circuit that generates them. In the previous section, the STN-GPe network was shown 239 to generate weak beta-band oscillations in the absence of exogenous beta inputs 240 (Fig. 2, 3, 4), and to phase lock to cortical beta-band inputs which amplified oscillatory 241 activity (Fig. 6). A potential role of the pallido-striatal loop could be to amplify

Mechanism of phase locking 288
The mechanism of phase locking of STN cells is illustrated in Fig. 8. Pooled cortical 289 spike trains (Fig. 8.A-B, green) show how sparse cortical beta bursts (Fig. 6.B) result in 290 distributed synaptic inputs to individual STN neurons that are not tightly phase locked, 291 but have a combined firing rate that is modulated at the beta frequency ( Fig. 6.B).

292
While these exogenous cortical inputs had high spike timing variability, STN and GPe 293 spikes became highly structured and tightly locked to the beta oscillation through the 294 feedback inhibition mechanism. The cortical beta modulation is transmitted to the STN 295 and then to the GPe through their excitatory projections (see phase vectors in 296 Fig. 7.Dii). When the inhibitory feedback arrives back in STN this shuts down spiking 297 April 11, 2019 10/33 ( Fig. 8.A) and simultaneously primes the cell for the next period of increased cortical 298 excitation by de-inactivating Ca 2+ channels ( Fig. 8.C) and N a + channels. As the 299 cortical firing rate rises again, synaptic currents ( Fig. 8.B) combine with dendritic Ca 2+ 300 currents to overcome any lingering inhibition and cause the next wave of phase-locked 301 STN spikes. The striatal beta inputs further decreased spiking variability of GPe 302 neurons by narrowing their time window of firing through phasic inhibition (purple 303 phase vector in Fig. 7.Dii). range for increased levels of drive (Fig. 2, 6.E). Besides affecting population firing rates 334 and oscillation frequencies, the excitation-inhibition balance also strongly influenced the 335 firing pattern of STN neurons: for a low ratio of excitation to inhibition and sufficiently 336 strong inhibitory currents, STN neurons transitioned to a firing mode characterized by 337 low-frequency tight bursts (high intra-burst firing rate, Fig. 2-4). Low-frequency 338 bursting was periodic at 2-5 Hz but was not synchronized between cells. This shift in 339 firing pattern towards sparse, tight bursting is in correspondence with changes in 340 burst-related measures like intra-burst firing rate and sub-beta band power that are 341 most predictive of akinetic-bradykinetic symptoms in humans [11] and monkeys [10], 342 although sub-beta band occurs at a higher frequency in monkeys. The firing rate and 343 pattern of GPe neurons was less sensitive than STN neurons to variations in its 344 excitatory or inhibitory drive as it was under negative feedback control by homeostatic 345 mechanisms that operated in synergy to stabilize its E/I ratio. However, GPe neurons 346 did synchronize more strongly under conditions of low excitatory drive from the STN 347 enabling them to act more autonomously and synchronize through inhibitory collaterals. 348 When beta-band spiking inputs were applied via cortico-STN afferents, the 349 STN-GPe network phase locked to the beta rhythm and input frequencies near the 350 autonomous oscillation frequency for a E/I ratio were preferentially amplified, reflected 351 in increased phase locking and oscillatory power of the somatic membrane voltage 352 ( Fig. 6.C,E). This is supportive of experimental observations that oscillatory activity in 353 STN is contingent on cortical oscillations [38], likely transmitted though the hyperdirect 354 pathway [12]. Phase locking and oscillation power were further strengthened by addition 355 of striatal oscillatory inputs with a particular phase relationship to cortical oscillatory 356 inputs (Fig. 7). The maximum in phase locking occurred when cortical and GPe promotes phase locking of STN neurons to beta-band cortical inputs [39].

369
Relation to other models of oscillations 370 The mechanism by which activity in the STN and GPe oscillates by means of 371 alternating phases of excitation and inhibition in a delayed negative feedback loop has 372 been described in previous models [19,21,48]. The mechanisms of oscillation in our 373 model of the STN-GPe network is consistent with this, but further illustrates the dual 374 role of GPe inhibition of STN neurons in shutting down ongoing activity and priming 375 for the next wave of excitation (Fig. 8). Moreover, in our model endogenously generated 376 oscillations were weak in the autonomous STN-GPe network in all simulation except 377 when fast-acting inhibitory currents in STN dominated strongly over slow ones. Spiking 378 variability was sufficiently decreased to cause a strongly synchronized beta-band 379 oscillation to emerge only when exogenous beta-range inputs were applied. (Fig. 6, 7). 380 April 11, 2019 12/33 The role of such patterning input could be fulfilled either by an extrinsic oscillator or by 381 reverberation of oscillations in connected feedback loops, e.g. pallido-striatal 382 circuits [24], recurrent striatal loops [23], or the thalamocortical 383 backprojection [25,28,56]. Our results, therefore, support the role of resonance with 384 other basal ganglia loops in generating exaggerated synchrony. Our model showed that 385 the preferred oscillation frequency could be modulated by varying the ratio of excitation 386 to inhibition in STN and GPe (Fig. 2B, 6.E) and that the ranges of strong oscillations 387 and low-frequency bursting overlapped (Fig. 6.C). Mean field models have shown that 388 the oscillation frequency is strongly modulated by transmission delays and membrane 389 time constants [21,57], while synaptic strengths have a weaker effect on the oscillation 390 frequency with sensitivity comparable to that in our model [21,22]. However, in a model 391 incorporating active ion channels that contribute to synaptic integration, synaptic ). In our model, low-frequency 400 fluctuations in beta-band power could arise from calcium-activated potassium currents 401 triggered by phases of high activity due to short-term plasticity dynamics as shown in 402 other models [30]. Under the full-loop resonance hypothesis low-frequency 'waxing and 403 waning' of beta activity could arise from modulations in beta power in any node of the 404 BGTC network, either through intrinsic dynamics or competition between extrinsic 405 inputs.

407
One of the main advantages of using a biophysically detailed model that incorporates 408 the level of detail presented here is that afferents from each pre-synaptic population can 409 target a specific region of the dendritic tree (Table 3, [61,62]. Hence, a model that incorporates a full 414 complement of ion channel and the synapse groups that interact with them may be 415 expected to yield a more realistic picture of how synchronization arises in the network. 416 In future studies, this can also contribute to a better understanding of neuronal currents 417 contributing to the local field potential in synchronized and asynchronous states, as 418 synaptic and ionic transmembrane currents combine to form the extracellular currents 419 that underpin this signal [63].

420
A second advantage of detailed models is that parameters have a clear relationship 421 to the underlying biophysical system and are more meaningful compared to the case 422 where parameters are lumped, as in single-compartment conductance-based models, or 423 abstracted as in mean-field or generalized integrate-and-fire models. This allows for 424 direct translation of experimental findings to parameter variations in the model. On the 425 other hand, detailed cell models are more sensitive to correct estimation of these 426 parameters which is limited by measurements performed for the purpose of model fitting 427 as well as the fitting procedure itself. In this model, downregulation of HCN channel 428 currents under dopamine depletion was modeled as a decrease in its peak conductance. 429 However, dopamine is known to interact with several ion channels that are involved in 430 linearizing the current-firing rate curve and regularizing autonomous pacemaking of 431 STN neurons [64][65][66] which are not all included in the STN cell model used here [43]. 432 Recent evidence suggests that this loss of autonomous spiking is a necessary condition 433 for the exaggerated cortical patterning of STN that is related to motor dysfunction [67]. 434 Better characterization of the ion channels involved in pacemaking and their response to 435 dopamine depletion will enable the systematic exploration of their contribution to STN 436 response properties and pathological firing patterns. based on their molecular profile and axonal connectivity [41]. Of these we only modeled 444 the prototypic population, projecting back to STN and preferentially firing in anti-phase 445 to it. Moreover, the GPe cell model used was only one representative candidate out of a 446 large set of models with varying ion channel expression and morphology that matched a 447 corresponding database of electrophysiological recordings [50]. Likewise, the STN model 448 represents a stereotypic characterization and does not capture variability in firing 449 properties and receptor expression. In particular, STN neurons in vivo are known to 450 have variable expression of GABA B receptors [44] which cause strong hyperpolarization 451 responses and longer pauses in some but not all STN neurons [52] and a strong rebound 452 burst response [44] in a subset of these. A model that takes into account the biological 453 variability in GABA B expression and that of channels underlying the rebound response 454 may reveal a wider range of responses to increased inhibition among STN neurons. In 455 such a model, beta rhythms could be transmitted to a subset of STN neurons whereas 456 others would show longer pauses with stronger rebound bursts. Moreover, the GABA B 457 synapse model used does not fully account for activation of extrasynaptic GABA B R due 458 to GABA spillover [44] which is mediated by tonic high-frequency and coincident firing 459 of afferents [40]. A model where multiple GABAergic synapses act on a shared pool of 460 extrasynaptic GABA B R might increase the importance of synchronized pre-synaptic 461 activity in switching STN neurons to a burst-firing mode.   Conductance-based models 512 In the conductance-based formalism, the membrane potential v j (mV) in each 513 compartment j of a multi-compartmental cable model is governed by: where x (cm) is the position along the cable, c m (µF/cm 2 ) is the specific membrane 515 capacitance, d (cm) is the cable diameter, R a is the specific axial resistance (Ωcm), g m 516 (S/cm 2 ) is the passive membrane conductance, E m (mV) the leakage reversal potential, 517 I ion,j (mA/cm 2 ) are the ionic currents flowing across the membrane of compartment j, 518 and I syn,j (mA/cm 2 ) are the synaptic currents at synapses placed in the compartment. 519 Each ionic current is governed by an equation of the form: where g x is the maximum conductance of the channel (S/cm 2 ), E x is the reversal with m ∞ (v) and τ m (v) representing the voltage-dependent steady state value and 524 time constant of the gate. For some currents the gating dynamics are described in terms 525 of the opening and closing rates α m and β m related through τ m = 1 αm+βm , 526 m ∞ = αm αm+βm : Reversal potentials are assumed constant unless otherwise noted. The reversal 528 potential for Ca 2+ currents was calculated using the Nernst equation from the intra-529 and extracellular ion concentrations: where T is the temperature in Kelvin, R is the universal gas constant, F is the 531 Faraday constant, and z is the valence of the Ca ion (+2). Intracellular calcium 532 buffering in a sub-membrane shell is modeled as: where c is a unit conversion constant, d is the thickness of the sub-membrane shell, 534 and τ Ca is the time constant of decay.

535
Synapses were modeled by a dual exponential profile with rise and decay times τ rise 536 and τ decay modulated by the fraction of synaptic resources in the active state which was 537 governed by Tsodyks-Markram dynamics [71]: where, g syn is the peak synaptic conductance, B-A represents the synaptic gating 539 variable, f peak is a normalization factor so that B-A reaches its maximum at time t peak 540 after the time of spike arrival t spk , R is the fraction of vesicles available for release, U SE 541 is the release probability, and τ rec and τ f acil are the time constants for recovery from The metabotropoc GABA B receptor-mediated current was modeled as an 547 intracellular signaling cascade based on the model by Destexhe and Sejnowski [73]. The 548 equations describing G-protein activation and the synaptic current were retained, but  Table 2. The channel density distributions are described extensively in ref. [43]. As a 566 source of noise, a current that followed a Gaussian amplitude distribution with mean 567 zero and standard deviation 0.1 was added to the somatic compartment.

568
The synaptic currents included an excitatory glutamergic input from cortex, acting 569 through AMPA and NMDA receptors, and an inhibitory GABAergic input from the 570 GPe, acting through GABA A and GABA B receptors (Table 3) NMDA receptors with majority NR2B and NR2D subunits that have dendritic punctual 579 expression [15]. The proximal synapses had only an NMDA component and represented 580 NMDA receptors with NR2A subunits that have faster kinetics [15]. The parameter 581 values for the conductance-based synapses are listed in Table 3.  Table 4. The channel density distributions are those described in ref. [50] for  The mean firing rate of STN surrogate spike sources was increased from 14.6 to 29.5 614 Hz in the parkinsonian state [2]. The peak GABA A and GABA B conductance of GPe to 615 STN synapses was increased by 50% and the GABA B decay time constant increased by 616 2 ms to model the increase in the number of contacts, vesicle release probability, and 617 decay kinetics of GPe afferents [74]. To model the reduction in cortico-STN axon 618 terminals and their dendritic targets [46,47] the number of CTX-STN afferents was 619 reduced to 70% of the normal condition, corresponding to the ratio of vGluT1 expression 620 in the normal and dopamine depleted condition used to label axon terminals [46]. To 621 model functional strengthening of remaining synapses, the AMPA and NMDA peak 622 conductances of remaining synapses were multiplied by the ratio of the current scaling 623 factors reported in [75] to the fraction of remaining synapses. Finally, HCN currents 624 were reduced by 50% to model reduced depolarization and spontaneous activity after 625 dopamine depletion [76,77,77] and modulation of HCN current by D2R receptors [65]. 626 In GPe neurons the peak AMPA conductance of STN afferents was increased by 50% 627 to model the modulatory effect of dopamine on glutamergic excitatory currents [78][79][80]. 628 The strengthening of GPe-GPe collaterals [81,82] was modeled by increasing the peak 629 GABA A and GABA B conductances by 50%. The mean firing rate of GPe surrogate 630 spike sources was decreased from 33.7 to 14.6 [9]. Finally, the HCN channel 631 conductance was decreased by 50% in accordance with experimental data [83].

632
Cortical projection neurons were modeled as Poisson spike generators firing at 10 Hz, 633 a multiple of the experimentally reported rate of 2.5 Hz [84], so that each synapse 634 represented the combined inputs of four pre-synaptic neurons (making use of the 635 additive property of the Poisson distribution). Where cortical oscillatory bursting 636 inputs were modeled, in each period of the beta oscillation 10% of cortical neurons were 637 selected randomly to emit a burst synchronously, using the same phase for the burst 638 onset but a variable number of spikes with inter-spike intervals sampled from the 639 interval [5,6] ms.