A biophysical model of striatal microcircuits suggests θ-rhythmically interleaved γ and β oscillations mediate periodicity in motor control

Striatal oscillatory activity is associated with movement, reward, and decision-making, and observed in several interacting frequency bands. Local field potential (LFP) recordings in rodent striatum show dopamine- (DA-) and reward-dependent transitions between two states: a “spontaneous” state involving β (~15-30 Hz) and low γ (~40-60 Hz), and a state involving θ (~4-8 Hz) and high γ (~60-100 Hz) in response to DAergic agonism and reward. The mechanisms underlying these rhythmic dynamics, their interactions, and their functional consequences are not well understood. In this paper, we propose a biophysical model of striatal microcircuits that comprehensively describes the generation and interaction of these rhythms, as well as their modulation by DA. Building on previous modeling and experimental work suggesting that striatal projection neurons (SPNs) are capable of generating β oscillations, we show that networks of striatal fast-spiking interneurons (FSIs) are capable of generating θ and γ rhythms. Our model consists of three interconnected populations of single or double compartment Hodgkin-Huxley neurons: a feedforward network of FSIs exhibits a D-type potassium current as well as DA-modulated gap junctional and inhibitory connectivity, and two networks of SPNs exhibit an M-type potassium current and express either excitatory D1 or inhibitory D2 DA receptors. Under simulated low DAergic tone the FSI network produces low γ band oscillations, while under high DAergic tone the FSI network produces high γ band activity nested within a θ oscillation. SPN networks produce β rhythms in both conditions, but under high DAergic tone, this β oscillation is interrupted by θ-periodic bursts of γ-frequency FSI inhibition. Thus, in the high DA state, packets of FSI γ and SPN β alternate at a θ timescale. In addition to a mechanistic explanation for previously observed rhythmic interactions and transitions, our model suggests a hypothesis as to how the relationship between DA and rhythmicity impacts motor function. We hypothesize that high DA-induced periodic FSI γ-rhythmic inhibition enables switching between β-rhythmic SPN cell assemblies representing the currently active motor program, and thus that DA facilitates movement by allowing for rapid, periodic shifts in motor program execution. Author summary Striatal oscillatory activity is associated with movement, reward, and decision-making, and observed in several interacting frequency bands. The mechanisms underlying these rhythmic dynamics, their interactions, and their functional consequences are not well understood. In this paper, we propose a biophysical model of striatal microcircuits that comprehensively describes the generation and interaction of striatal rhythms, as well as their modulation by DA. Our model suggests a hypothesis as to how the relationship between DA and rhythmicity impacts the function of the motor system, enabling rapid, periodic shifts in motor program execution.

Behavior of single model FSI over a range of inputs and D-current conductances. (A) i. A single model FSI with low tonic excitation (8µA/cm 2 ) spikes at a low γ frequency nested in slow bursting, while a single model FSI with high tonic excitation (20µA/cm 2 ) spikes at a high γ nested in slow bursting. ii. Power spectral density of voltage traces in (A)i, comparing low and high levels of tonic excitation. (B) i. Single model FSI with tonic excitation (10µA/cm 2 ) and weak Poisson noise (λ = 500) spikes at γ nested in θ, while a single model FSI with tonic excitation (10µA/cm 2 ) and strong Poisson noise (λ = 5500) has limited low-frequency content. ii. Power spectral density of voltage traces in (B)i, comparing low and high levels of noise. (C) Three-dimensional false-color plot demonstrating the dependence of the bursting regime on g d and I app .
(D) Three-dimensional false-color plot demonstrating the dependence of firing rate on g d and I app .
quiescence to (periodic) bursting to periodic spiking. The bursting regime, of particular 69 interest in this work, is dependent on the level of tonic excitation (Fig. 1A,C,D), the 70 level of noise in excitatory drive ( Fig. 1B & 2C), and, centrally, the D-current 71 conductance (Fig. 1C,D). With lower levels of D-current (as used in previous FSI 72 models [10,17,19]), bursting is aperiodic. For sufficiently large D-current conductance, 73 FSI bursting occurs for a broad range of applied currents (I app over 8 uA/cm2, Fig.   74 1C,D). The frequency of bursting depends on the decay time constant of the D-type 76 potassium current (τ D ); in the absence of noise, it is in the δ frequency range for 77 physiologically relevant τ D (<∼200 ms, Figure 1A,C, 2B). Note that τ D changes the 78 inter-burst interval without changing the timing of spikes within a burst (Fig. 2B). 79 Burst frequency also increases with increasing tonic excitation (Fig. 1,A,C), and 80 increases to θ frequencies with small amounts of noise (Fig. 1B, 2C), which decrease the 81 interburst interval. However, large amounts of noise abolish rhythmic bursting 82 altogether (at least in single cells, Fig. 1B). 83 As shown previously [17], the FSI model's γ rhythmic intraburst spiking arises from 84 its minimum firing rate, which is also set by the D-current conductance ( Fig. 2A); when 85 this conductance is zero, the model has no minimum firing rate ( Fig. 2A). As this conductance is increased, the minimum firing rate also increases. Thus, our choice of 87 D-current reflects not only our interest in the bursting regime, but also our desire to 88 match experimental observations of striatal γ frequency [2,17]. FSI spiking frequency 89 also increases with tonic drive ( Figure 1A,D). Since simulated DA acts on our FSI 90 model by increasing tonic excitation, DA causes a switch in model FSI spiking from low 91 γ rhythmicity to δ-modulated high γ rhythmicity. Below, we demonstrate that the FSI 92 γ is determined by this single-cell rhythmicity and is mostly independent of the 93 timescale of inhibitory synapses. 94 In summary, a single model FSI displays low-frequency-nested γ oscillations, 95 dependent on the D-type current, under a wide range of tonic excitation levels. Both low 96 frequency power and γ frequency increase with tonic excitation. While noise increases 97 the frequency of the slower rhythm from δ to θ, it also diminishes the power of this 98 rhythm in the single cell. Below we demonstrate that all of these effects are also present 99 in a network of FSIs, with the key difference that the θ rhythm becomes robust to noise. 100 FSI networks produce DA-dependent θ and γ rhythms 101 To determine if θ and γ oscillations persist in networks of connected FSIs, and how DA 102 could modulate network dynamics, we simulated a network of 100 model FSIs connected 103 randomly (with an independent connection probability of 0.3 for each type of 104 connection) by both inhibitory synapses and gap junctions. We also implemented three 105 salient and experimentally observed effects of DA on FSI networks: increased tonic 106 excitation of individual FSIs [12], increased gap junction conductance between FSIs [13], 107 and decreased inhibitory conductance between FSIs [12] (see Methods).

108
Unlike in single cells, FSI network θ rhythmicity is dependent on sufficient levels of 109 tonic excitation: at low levels of tonic input (I app <≈ 1µA/cm 2 ), the FSIs do not attain 110 enough synchrony for a strong network θ (Fig. 3Aii). As in single cells, FSI network θ 111 power increases with tonic input strength (Fig. 3Aii). Sufficiently strong gap junction 112 coupling is also a requirement for the FSI network to attain sufficient synchrony to 113 produce θ rhythmicity (Fig. 3Cii), protecting the FSI network θ rhythm from the effects 114 of noise (as in [14]). Finally, inhibitory synaptic interactions between FSIs have a 115 desynchronizing effect that interferes with network θ, and increasing inhibitory conductance within the FSI network decreases power in the θ band (Fig. 3Bii). Power and frequency of θ and γ rhythms in FSI network mean voltage as a function of (A)tonic input current, (B) GABA A conductance, and (C) gap junction conductance.. The parameters not being varied in plots A-C are held at the high DA values (g gap = 0.2mS, g syn = 0.005mS, I app = 16uA/cm2, tau gaba = 13ms . (D) Gamma frequency as a function of GABA a synaptic time constant and level of dopamine.
FSI network γ power and frequency both increase with tonic input strength ( Fig.   118 3Ai), and, like the network θ, the network γ rhythm is dependent on sufficient gap 119 junction conductance and is disrupted by inhibition ( Fig. 3B & C, i).

120
To explore FSI network dynamics that might be observed during normal fluctuations 121 in DA during goal-directed tasks [20], we simulated FSI network activity under two 122 conditions, simulated low (or baseline) and high DAergic tone. During simulated low 123 DAergic tone, characterized by low levels of FSI tonic excitation and gap junction 124 conductance, and high levels of inhibitory conductance, the network produces a 125 persistent low frequency γ oscillation (∼ 55 Hz) in the model LFP (the mean voltage of 126 the FSI network, Fig. 4Bi-Di). The raster plot of FSI spike times (Fig. 4Eii) shows that 127 individual FSIs exhibit sparse spiking in the low DA state. Although individual FSIs 128 exhibit periodic spike doublets or bursts (γ-paced and entrained to the network γ) that 129 recur at θ frequency, the timing of these bursts is independent (Fig. 4Ei). Therefore, 130 while θ power is present at the level of individual FSIs, there is not sufficient synchrony 131 for it to appear in the network (Fig. 4Di).    SPNs were connected by all to all inhibitory synapses within (connection probability 1) 176 but not across populations. There were no connections from SPNs back to FSIs [22]. Prior work has shown γ oscillations in striatal FSIs arising from an interaction between 206 the spiking currents and the spike frequency adaptation caused by the potassium 207 D-current, which produces a minimum FSI firing rate in the γ range [17,24]. The 208 frequency of the FSI γ depends on excitatory drive to the FSIs, which in our model 209 leads to the modulation of γ frequency by DA, a phenomenon also observed in striatal γ 210 oscillations in vivo [25][26][27][28].

211
Prior work has also suggested that the D-current is responsible for the bursting or 212 stuttering behavior of FSIs, in which brief periods of high frequency activity are 213 interspersed with periods of quiescence [10]. However, regularities in these periods of 214 quiescence have not been previously observed. Thus, the present study is novel in its 215 description of the generation of low-frequency rhythms by FSIs with high levels of D-current conductance; FSIs have previously been characterized solely as generators of 217 γ oscillations. Our model predicts that FSI-mediated slow rhythms depend on a high 218 level of D-current conductance. In our model, the D-current is activated by burst 219 spiking, e.g., at γ frequency, and hyperpolarizes the cell for roughly a θ period due to its 220 long time constant of inactivation. Though individual cells produce a δ rhythm, the 221 frequency of the resulting θ oscillation in the network is robust to changes in excitatory 222 drive. This transition to a higher rhythm in the network is likely a result of 223 gap-junction induced synchrony driving burst frequency higher while maintaining 224 robustness to noise. Notably, this study is also a novel demonstration of the generation 225 of both θ and γ oscillations by a single membrane current. condition is bursting sufficiently synchronized that θ power is present in the network.

233
The presence of θ at the network level can be attributed to the higher level of gap 234 junction conductance in the high DA condition (Fig. 3Cii). 235 The ability of gap junctions to generate synchrony is well established in 236 computational work [29]. Previous models from other groups suggest that gap junctions 237 can enable synchronous bursting in interneurons, such that the burst envelopes are 238 aligned, as in our model [14]. While a shunting effect of low conductance gap junctions 239 can inhibit spiking [30], gap junctions with high enough conductances have an 240 excitatory effect, promoting network synchrony [31,32]. Previous work has also shown 241 the importance of gap junction connectivity in stabilizing network γ oscillations in 242 silico [33], as well as network γ and θ oscillations in inhibitory networks in vitro and in 243 silico containing noise or heterogeneity [32]. FSIs in vivo are highly connected by gap 244 junctions as well as inhibitory synapses, [34] similar to the networks of inhibitory 245 interneurons that produce ING rhythms [35]. Unlike ING, however, our FSI network γ 246 is independent of GABAergic synapses: inhibitory conductance has only a small impact 247 on γ frequency, and γ power is highest when inhibitory synapses are removed ( Figure   248 3B). In slice, the γ resonance of striatal FSIs is dependent on gap junctions but not on 249 GABA [36], suggesting that our model is an accurate representation of FSI γ.

250
It is important to note that while our model is conceived as a representation of the 251 dorsal striatal circuit, physiologically similar fast-spiking interneuron networks are 252 present in cortex [10]; therefore, the mechanisms described here may contribute to the 253 generation of θ-modulated γ oscillations in cortex as well. ventral striatum in vivo (θ, β, low γ, and high γ) [37]. Previous modeling and 257 experiments suggest β can be generated by striatal SPNs [15,38,39]. Our results suggest 258 that FSIs generate striatal γ, and that motor-and reward-related increases in γ power 259 reflect increased striatal FSI activity.

260
There is evidence to support the existence of a locally generated striatal γ oscillation 261 that is not volume conducted and that responds to local DAergic tone [40,41]. The FSIs 262 of the striatum are the most likely candidate generator of this rhythm: they are unique 263 among striatal cell types in preferentially entraining to periodic input (from each other 264 and from cortex) at γ frequencies [42][43][44][45][46]. Different populations of striatal FSIs in vivo 265 entrain to different γ frequencies, and FSIs entrained to higher frequencies are also more 266 entrained to cortical input [25][26][27][28]. It is likely that different subpopulations of FSIs 267 selectively entrain to specific γ frequencies depending on physiological differences, 268 context, and neuromodulatory (e.g. DAergic) states; the frequency of γ may itself 269 determine cell assembly size and membership [33]. suggesting its responsiveness to DA (known to phasically increase in response to reward 275 cues) [4,[52][53][54]. θ has also been shown to modulate the response of SPNs to reward [55]. with both [1,[3][4][5][6]56]. β oscillations in the basal ganglia are thought to provide a "stay" 280 or "status quo" signal that supports maintenance of the currently active motor 281 program [16], and they are causally implicated in motor slowing and cessation [8,[56][57][58][59][60] or "program on" state marked by SPN β oscillations, and a "switch" or "program off" 292 state marked by FSI high γ oscillations, and that the θ period limits the speed of 293 sequential motor program execution (Fig. 6E).

294
In support of this hypothesis, striatal representations of behavioral "syllables" that 295 can be combined to create motor programs are active for a maximum of ∼200 ms [61], 296 and the velocity of continuous motion is modulated intermittently at a θ frequency 297 (∼6-9 Hz) [62]. In healthy animals, the duration of β bursts has an upper limit of ∼120 298 ms, about half a θ cycle [8], in agreement with our hypothesis of θ phase-modulation of 299 β activity. Striatal γ has also been observed in transient (∼150 ms) bursts that are 300 associated with the initiation and vigor of movement [63]. Additionally, other 301 biophysically constrained computational models have suggested that SPN assemblies 302 fire in sequential coherent episodes for durations of several hundred milliseconds, on the 303 timescale of one or several θ cycles [64]. Overall, evidence supports the hypothesis that 304 β and γ oscillations in striatum in vivo, and therefore the motor states they encode, are 305 activated on θ-periodic timescales.

306
Furthermore, β and γ power are anticorrelated in EEG and corticostriatal LFP [6,7,65], in agreement with our model's prediction that these rhythms are coupled 308 to opposite phases of ongoing θ rhythms. FSI and SPN firing are inversely correlated in 309 vivo, entrained to θ, and they are active during opposite phases of θ, as observed in our 310 model [23,48,[66][67][68]. θ-γ cross-frequency coupling is observed in striatum and increases 311 during reward, when DAergic tone is expected to be high [7,[69][70][71][72]. Our model suggests 312 that these cross-frequency relationships occur in part due to FSI inhibition of SPNs.

313
Though FSIs are smaller in number, FSI-SPN synapses have a much stronger effect 314 than SPN-SPN connections, with each FSI inhibiting many SPNs [22,73]. GPe [19], also based on an earlier model of stuttering FSIs [10]. In contrast to this work, 333 we emphasize the mechanisms producing β and the coordination of β and γ by θ, not 334 addressed previously [19].

Implications for disease
In Parkinson's disease, which is characterized by motor deficits and chronic DA depletion, β power is correlated with the severity of bradykinesia [1]. Parkinsonian β 338 may be generated by striatal D2 SPNs [15,38,39]. Parkinsonian conditions also produce 339 high cholinergic tone [75], known to decrease the conductance of GABAergic FSI-SPN 340 synapses [76]. Thus, the failure of the FSI inhibition-mediated motor program switching 341 described above may play a role in the motor deficits observed in Parkinson's: if DA is 342 low, and FSIs are unable to inhibit either D1 or D2 SPNs, θ modulation of SPN β 343 rhythmicity will be supplanted by ongoing D2 β rhythmicity, impairing motor initiation 344 by reducing the possibility of motor program switching in the Parkinsonian striatum.

345
In hyperkinetic motor disorders, γ and θ rhythms are potentiated: a mouse model of 346 Huntington's disease (HD) displays unusually high θ and γ band striatal LFP 347 power [77][78][79]; and L-DOPA-induced hyperkinetic dyskinesia is also characterized by 348 increased high γ and θ power and reduced β power in the striatal LFP [80][81][82][83] program switching [89]. A reduction in theta/γ cross frequency coupling has been 361 reported in L-DOPA-induced dyskinesia, suggesting that a chronic hyperkinetic high-DA 362 state may also abolish the FSI-generated θ-coupled γ produced here, possibly by 363 pushing the FSI out of its bursting regieme and into a tonic spiking mode [90]. These

381
In general, many DA-dependent changes in striatal neurophysiology have been 382 observed. For the sake of simplicity, most of these have been left out of our modeling.

383
For example, D1 and D2 SPNs respond differently to adenosine [91] and peptide 384 release [92], but we did not consider these significant factors in the production of 385 striatal β oscillations. 386 We also omitted inhibitory connections between D1 and D2 SPN populations. The 387 connectivity from D1 to D2 SPNs is very sparse (6 percent in vivo, the participation of multiple FSI populations is likely coordinated by cortex.

394
Coordinated FSI activity has proven hard to observe over long periods in vivo [28,93]. 395 However, FSIs form local functional circuits [94], and in vivo, striatal FSI assemblies as opposed to the overlapping projections modeled in our current study, and these 400 distinct populations respond differently to cortical oscillations [49]. Thus, local γ 401 synchrony may exist in small striatal subnetworks and be amplified by DA or cortical 402 input via the recruitment of multiple FSI subpopulations.

403
Finally, cortical input to both FSIs and SPNs was simulated as Poisson noise. In a 404 sense, then, we simulated a model of striatum to which cortex is not providing 405 informative input. It could be the case that this is a population that is not "selected" 406 by cortex to take part in motor activity, a population that is in a "listening" state 407 awaiting cortical input, or a population taking part in a learned behavior that can be 408 executed without cortical input. However, cortical input is probably essential in All membrane currents have Hodgkin-Huxley-type conductances formulated as: Each current in Eqn.2 has a constant maximal conductance (ḡ) and a constant (I m ) [11]. We do not model SPN up and down states which are not prevalent in the 473 awake state of striatum [96], the state being modeled, and therefore we do not include 474 the Kir current in our model, which is active during the SPN down state.

475
The sum of all excitatory inputs from the cortex and thalamus and inhibitory inputs 476 from striatal interneurons is introduced into the model using a background excitation 477 term (I app ). I app is the sum of a constant term and a Gaussian noise term.
The maximal fast potassium channel conductance isḡ K = 80mS. The reversal 498 potential for potassium is E K = −100mV .  2011 [15] and is the only synaptic connection between SPNs and from FSIs to SPNs.

520
The GABA A current has a single activation gate dependent on the pre-synaptic voltage. 521 The maximal GABA A conductance between FSIs isḡ ii = 0.08mS. The maximal 522 GABA A conductance from FSIs to SPNs isḡ ii = 0.6mS and between SPNs was 523 g ii = 0.1mS. These values are consistent with FSI to SPN inhibition being 524 approximately six times stronger than inhibition between SPNs [9].

525
The gating variable for inhibitory GABA A synaptic transmission is represented by s i . 526 For the j th neuron (FSI or SPN) in the network: The variable S i k ij describes the kinetics of the gating variable from the k th 528 pre-synaptic neuron to the j th post-synaptic neuron. This variable evolves in time 529 according to: The GABA A time constant of decay (τ i ) is set to 13 ms for SPN to SPN 531 connections [97] as well as for FSI to FSI connections and FSI to SPN connections [30] 532 The GABA A current reversal potential (E i ) for both FSIs and SPNs is set to -80 mV. 533 The rate functions for the open state of the GABA A receptor (g GABAA (V k )) for SPN to 534 SPN transmission is described by: The rate functions for the open state of the GABA A receptor (g GABAA (V k )) for FSI 536 to FSI and FSI to SPN transmission is: The value of τ r is 0.25 ms. FSIs were additionally connected by dendritic electrical 538 connections. The electrical coupling for dendritic compartment i is denoted as I elec , has 539 units in µA/cm 2 and is formulated as: The value of the gap junction conductance g gap depended on DA level (see below). 541 Within the 100-cell FSI network, each pair of FSIs had an independent 30 percent 542 chance of a dendro-dendritic gap junction chosen from a uniform random distribution 543 and an independent 30 percent chance of a somato-somatic inhibitory synapse also 544 chosen from a uniform distribution. SPNs are connected with each other in a mutually 545 inhibitory GABAergic network [100]. We modeled all to all connectivity of inhibitory 546 synapses from any SPN to any SPN of the same receptor subtype.

547
DA 548 DA impacts both connectivity and excitability in the model networks. DAergic tone was 549 simulated as having five components: direct excitation of FSIs [12], increased gap junction conductance between FSIs [13], decreased inhibitory conductance between FSIs [12], increased excitation to D1 SPNs, and decreased excitation to D2 SPNs.

552
DA-induced changes to SPN excitation were discussed above. Excitation to FSIs was 553 modeled as the sum of a tonic, DA dependent input current (I tonic ) and a noise term.

554
DA did not change the noise term in either SPNs or FSIs. The baseline DAergic tone 555 state was modeled in FSIs using I tonic = 4µA/cm 2 , g gap = 0.05mS and the GABA A 556 conductance between FSIs was g ii = 0.1mS. The high DA state was modeled in FSIs 557 using I tonic = 14µA/cm 2 , g gap = 0.2mS and g ii = 0.005mS.

558
Local field potential 559 The local field potential (LFP) was calculated as the sum of all voltages in all cells.

560
Stationarity of the network appears in the raster plots after about 500 ms. To eliminate 561 transients due to initial conditions, our LFP is evaluated only after 1,000 ms of 562 simulated time. We estimated the power spectral density of the simulated LFP using 563 the multitaper method [101].

564
Simulations 565 All simulations were run on the MATLAB-based programming platform DynaSim, a 566 framework for efficiently developing, running and analyzing large systems of coupled 567 ordinary differential equations, and evaluating their dynamics over large regions of 568 parameter space [102]. DynaSim is open-source and all models have been made publicly 569 available using this platform. All differential equations were integrated using a 570 fourth-order Runge-Kutta algorithm with time step was .01 ms. Plotting and analysis 571 were performed with inbuilt and custom MATLAB code.

572
Acknowledgments 573 We thank Jason Sherfey and Erik Roberts for developing the DynaSim Toolbox used to 574 run the simulations in this manuscript and assisting with debugging and visualization. 575 In addition, we thank Sean Patrick for assistance in interpreting the behavioral 576 implications of our model.