Abstract
The cortico-basal ganglia circuit is needed to suppress prepotent actions and to facilitate controlled behavior. Under conditions of response conflict, the frontal cortex and subthalamic nucleus (STN) exhibit increased spiking and theta band power, which are linked to adaptive regulation of behavioral output. The electrophysiological mechanisms underlying these neural signatures of impulse control remain poorly understood. To address this lacuna, we constructed a novel large-scale, biophysically principled model of the subthalamopallidal (STN-globus pallidus externus) network and examined the mechanisms that modulate theta power and spiking in response to cortical input. Simulations confirmed that theta power does not emerge from intrinsic network dynamics but is robustly elicited in response to cortical input as burst events representing action selection dynamics. Rhythmic burst events of multiple cortical populations, representing a state of conflict where cortical motor plans vacillate in the theta range, led to prolonged STN theta and increased spiking, consistent with empirical literature. Notably, theta band signaling required NMDA, but not AMPA, currents, which were in turn related to a triphasic STN response characterized by spiking, silence, and bursting periods. Finally, theta band resonance was also strongly modulated by architectural connectivity, with maximal theta arising when multiple cortical populations project to individual STN “conflict detector” units because of an NMDA-dependent supralinear response. Our results provide insights into the biophysical principles and architectural constraints that give rise to STN dynamics during response conflict, and how their disruption can lead to impulsivity and compulsivity.
SIGNIFICANCE STATEMENT The subthalamic nucleus exhibits theta band power modulation related to cognitive control over motor actions during conditions of response conflict. However, the mechanisms of such dynamics are not understood. Here we developed a novel biophysically detailed and data-constrained large-scale model of the subthalamopallidal network, and examined the impacts of cellular and network architectural properties that give rise to theta dynamics. Our investigations implicate an important role for NMDA receptors and cortico-subthalamic nucleus topographical connectivities in theta power modulation.
Introduction
Various studies have implicated the cortico-subthalamic “hyperdirect-pathway” in response inhibition on encountering response conflict (Baunez and Robbins, 1997; Baunez et al., 2001; Nambu et al., 2002; Frank, 2006; Aron et al., 2007; Isoda and Hikosaka, 2007, 2008; Verbruggen and Logan, 2008; Cavanagh et al., 2011; Zaghloul et al., 2012; Schmidt et al., 2013; Zavala et al., 2014; Herz et al., 2016; Kelley et al., 2018; Wessel et al., 2019). The STN mediates these functions via global suppression of motor output, putatively via diffuse projections to basal ganglia (BG) output structures (Parent and Hazrati, 1995; Frank, 2006; Wessel et al., 2019). Disruption of subthalamic nucleus (STN) causes impulsive behavior (Baunez et al., 1995, 2001; Baunez and Robbins, 1997; Frank et al., 2007; Wylie et al., 2010; Coulthard et al., 2012; Green et al., 2013; Fife et al., 2017; Ghahremani et al., 2018). Accounting for these observations, computational models of BG (Frank, 2006; Bogacz and Gurney, 2007; Ratcliff and Frank, 2012; Wiecki and Frank, 2013) have further posited that elevated STN firing dynamically raises a “decision-threshold” to afford time to resolve the conflict.
Despite the growing body of work on STN dynamics, the electrophysiological mechanisms remain poorly understood. Many studies have reported elevated theta (4-8 Hz) power in frontal cortex and STN during response conflict (Cavanagh et al., 2011, 2014; Zavala et al., 2013, 2014, 2016; Herz et al., 2016). Such theta modulations induce transient global suppression of motor activity (Wessel et al., 2019), leading to slower but more accurate choices, as captured by elevated decision thresholds in quantitative models (Cavanagh et al., 2011; Zavala et al., 2014, 2016; Herz et al., 2016, 2017; Kelley et al., 2018). But what are the mechanisms of such STN theta band signals, and how are they modulated? For example, while STN theta correlates positively with decision threshold and response time slowing during high-conflict conditions, the opposite has been observed during low conflict (Cavanagh et al., 2011; Herz et al., 2016), thereby questioning whether theta power per se indicates the need for response caution. Moreover, many factors could influence STN theta, from conductance dynamics of individual cells (Bevan and Wilson, 1999; Cooper and Stanford, 2000; Wilson, 2010; Deister et al., 2013; Jones, 2016; Rubin, 2017) to network topography (Mathai and Smith, 2011; Nambu, 2011; Haynes and Haber, 2013; Schroll and Hamker, 2013; Kita et al., 2014; Alkemade et al., 2015), to the dynamics of the cortical inputs to STN (Zavala et al., 2014).
To query the conditions under which STN theta power is modulated, we developed a novel biophysically principled large-scale model of the subthalamopallidal circuit by adapting previous STN and globus pallidus externus (GPe) cellular models (Terman et al., 2002; Rubin and Terman, 2004). Our model incorporates arkypallidal GPe units and imposes sparse and probabilistic connectivities (Oorschot, 1996; Bevan et al., 2007; Kita, 2007, 2010; Sadek et al., 2007; Baufreton et al., 2009; Wilson, 2010; Mallet et al., 2012; Abdi et al., 2015). We first constrained our network to capture single-cell electrophysiological patterns. We then characterized the effects of various dynamic regimens of cortical inputs (Jones, 2016) on STN theta, and how they are altered by cellular mechanisms (e.g., AMPA and NMDA currents). We simulated response conflict in terms of coactivation of multiple cortical populations (representing competing motor responses), and explored the impact of their topographical connectivities on STN subpopulations. Finally, we investigated how these mechanisms affect both oscillatory and spiking dynamics in these subpopulations.
Our simulations implicate NMDA-dependent mechanisms within STN in conflict-induced elevations in theta power and spiking. Theta band resonance was also strongly modulated by architectural constraints, with maximal response achieved when multiple cortical inputs converged on STN subpopulations. Analysis of the underlying mechanism revealed an NMDA-dependent supralinear response in STN “conflict-detector” units. These findings provide a potential resolution to various conundrums in the empirical literature and a mechanistic basis to guide potential therapeutic developments for impulsive behaviors.
Materials and Methods
Overview of subthalamopallidal (STN-GPe) network
We constructed a large-scale model of the STN-GPe network (see Fig. 1A) to study the conditions regulating STN theta and spiking during response conflict. We began by adapting prior models derived from rat electrophysiology of the STN and GPe (Terman et al., 2002; Rubin and Terman, 2004), expanding them to include the following: (1) two distinct subpopulations of GPe, (2) NMDA currents, and (3) various types of cortical inputs. STN activity is predominantly patterned by reciprocal excitatory and inhibitory interactions with the GPe (Baufreton et al., 2009; Wilson, 2010, 2013), and we thus considered it important to incorporate the recently established arkypallidal (GPeA) subclass dynamics and connectivities (Mallet et al., 2012; Abdi et al., 2015), in addition to the prototypic pallidal (GPeP) cells. Because we focused on the subthalamopallidal network, omitting the larger circuit (e.g., striatal input, globus pallidus internus output) in which it is embedded, we began by tuning the baseline state of the network (its intrinsic state with no cortical drives; see Baseline tuning of the STN-GPe network), to reflect in vivo spiking statistics. While capturing the overall effects of the nonmodeled nuclei, this trade-off facilitated our model tuning to be challenged by biophysical constraints (see Biophysical constraints and network size) of the STN and GPe units at various levels.
Our simulations were implemented using the Neuron framework in Python with Message Passing Interface Parallelization and Parallel Context (Hines et al., 2009; Jones et al., 2009), and were run on High Performance Computing resources at Brown University and XSEDE resources (Sivagnanam et al., 2013; Carnevale et al., 2014; Towns et al., 2014). Data were stored in the HDF5 (Folk et al., 2011) hierarchical structure with Parallel Input/Output functionality.
Biophysical constraints and network size
Spanning various levels of neuroscience, this study demanded that mechanisms regarding synaptic kinetics, cellular dynamics, intranuclei and internuclei connectivities conform to biophysical constraints before investigating how these mechanisms may give rise to neural signatures of response conflict. Given the inherent complex nonlinearity across and within several scales, and the compounding tight coupling of biophysical objectives, tuning such a network necessitated several iterative steps (see Fig. 1B, as described in Tuning synaptic dynamics).
Electrophysiological constraints
We first imposed that postsynaptic currents reproduce the rise and decay kinetics reported under voltage-clamp experiments, and that synaptic conductances show voltage dependence (see Fig. 1G). At the cellular level, we ensured that membrane potentials exhibit physiologically plausible dynamics related to action potential (AP) generation and resting membrane potentials (RMPs); those that looked pathologic (e.g., elevated RMP, APs with broad spike widths and hyperpolarizing peak voltages) were rejected. We also insisted that each subclass of cells accorded with in vivo population spiking statistics in terms of means and distribution shapes (see Fig. 1E) subject to connectivity constraints. Furthermore, because synchronous STN and GPe firing is characteristic of pathologic conditions (Goldberg et al., 2004; Hamani et al., 2004; Hammond et al., 2007; Eusebio et al., 2009; Wilson, 2010, 2013; Schwab et al., 2013), minimizing the incidence of such synchrony was an additional biophysical constraint on the target baseline state, which has also been captured by other computational models (Hahn and McIntyre, 2010; Park et al., 2011; Pavlides et al., 2015). These requirements, coupled with the need to have physiologically relevant temporal resolution to model theta signals, motivated our use of active conductance Hodgkin-Huxley based models (Terman et al., 2002; Rubin and Terman, 2004) (for somata, see Cellular tuning). Because available data from paired recordings were measured at somatic levels only, with no information as regards dendritic ion channel distributions and postsynaptic potential attenuations, we chose single compartment models.
Anatomical constraints
The sparsity in subthalamopallidal connectivities (Baufreton et al., 2009; Wilson, 2010) (see Fig. 1A) required a large network size. However, the stiff dynamics of the GPe cells dictated a very small simulation step size of 0.0025 ms for numerical accuracy, which placed an upper limit on the network size given our computational resources. These conflicting restrictions, in conjunction with reported anatomic population ratios of STN to GPe being one to three across several species (Oorschot, 1996; Hardman et al., 2002), informed our network size of 300 STN and 900 GPe cells. Moreover, as GPeP cells are estimated to be around twice the population size of GPeA cells (Mallet et al., 2012; Abdi et al., 2015), we divided the GPe population into 600 GPeP and 300 GPeA units.
Overview of the tuning process
Figure 1B provides a flowchart of the model tuning process. We started by tuning the synaptic kinetics to reproduce the rise and decay of postsynaptic currents under voltage-clamp experiments recorded at the somata (see Synaptic tuning). Using the GPe model in Rubin and Terman (2004) as a template, which would correspond to the prototypic subclass in our case, we changed the intrinsic properties to derive the GPeA model. Fixing the connection probabilities of the STN, GPeP (Baufreton et al., 2009; Wilson, 2010) and our new GPeA (as in Fig. 1A), we then performed a thorough search in the synaptic conductance space and evaluated whether the baseline states of the generated networks met the imposed in vivo spiking statistics (see Fig. 1E). Unfortunately, we found no parameter sets from the original Rubin and Terman (2004) model that could conform to the latter constraints. Consequently, we tuned the dynamics of all the cellular units, namely, STN, GPeP, and GPeA (i.e., departing from those used by Rubin and Terman, 2004), and repeated the synaptic conductance space search and evaluations. The latter paired tuning process of cellular dynamics followed by synaptic conductance was iterated until a parameter set that met several biophysical constraints was obtained. We note here that several parameter sets, while successful in capturing the desired spiking statistics, were eventually rejected when challenged by biophysical constraints at the cellular level, such as AP characteristics (see Cellular tuning).
Tuning synaptic dynamics
Lack of data to tune metabotropic postsynaptic influences on single cells restricted our investigations to ionotropic currents (Eq. 1a) only.
In contrast to solving the first-order differential model for synaptic activity in Rubin and Terman (2004) and Terman et al. (2002), we opted for the computationally less expensive biexponential model with rise and decay constants, Neuron's Exp2Syn mechanism (Carnevale and Hines, 2006), Equation 1b, for synaptic responses: we modified to allow maximal conductances
Parameters for synaptic modelsa
Postsynaptic current:
Closed-form solution for biexponential model of synaptic conductance:
Voltage-dependent maximal synaptic conductance model:
Despite its simplicity, the biexponential model does not have a closed-form solution to solve for the (10%-90%) rise τrise and decay τdecay rates. To replicate the time courses of synaptic responses under voltage clamps, we manually fitted the biexponential model by guessing an initial solution of time rises and decays to achieve the 10%-90% characteristics, iteratively refining the ranges to converge on a satisfying solution, which are reported in Table 1.
Individual cell dynamics and synaptic connections
STN and GPeP cells were modeled as in Terman et al. (2002), with the Na+, K+, leak,
All equations and parameters for STN and GPeP were as in Rubin and Terman (2004) and Terman et al. (2002), up to specific capacitance Cm scaling to 1 μF/cm2 for every cell, except for changes detailed in Table 2.
Parameters for cellular spontaneous activitya
Membrane dynamics:
AP generating currents:
Ca2± related currents:
Postsynaptic currents:
First-order kinetics for gating variables: n, h:
Steady-state voltage dependence for gating variables: m, n, and h:
Spike detection thresholds for activating the NetCons were set at −47.4, −56.6, and −55.0 mV for STN, GPeP, and GPeA units, respectively.
Each STN unit had an AMPA synapse and an NMDA synapse responding to cortical spike trains. The cortico-STN synaptic parameters (Table 1) were matched to data from Chu et al. (2015), both for kinetics and maximal conductance voltage dependence. The NMDA to AMPA maximal conductance ratio was kept at a fixed value corresponding to experimentally measured NMDA conductance at 40 mV divided by the AMPA conductance at −80 mV, which in this case is 1.402. We use fscale = 1 for NMDA at 40 mV and AMPA at −80 mV in the voltage-dependent maximal conductance formalism (Eq. 1e; see Fig. 1G); hence, the maximal conductances would match the experimental values. Each STN unit had a GABA synapse for responding to GPeP stimulation. Both GPeP and GPeA cells had an AMPA and an NMDA that received connections from the STN, and a GABA synapse that was triggered by GPeP and GPeA activities. The different units were probabilistically connected as shown in Figure 1A, with the values in square boxes representing the sparse connection probabilities, based on previous reports (Bevan et al., 2007; Kita, 2007, 2010; Sadek et al., 2007; Baufreton et al., 2009; Wilson, 2010; Mastro et al., 2014; Abdi et al., 2015). For every pair of presynaptic and postsynaptic cells, we performed a Bernoulli trial based on the predefined sparse probabilities, and connected the pair if the result was 1.
Baseline tuning of the STN-GPe network
After probabilistically connecting the network, we first tried to reproduce a baseline state that represents a healthy condition, where the units in the network fire asynchronously within and across subpopulations (see Fig. 1C–F), with normally distributed spiking rates with means for STN, GPeP, and GPeA set at ∼11 (Kreiss et al., 1997; Beurrier et al., 1999; Bevan and Wilson, 1999), 30 and 5 spikes/s respectively (Kita and Kitai, 1991; Nambu and Llinas, 1994; Cooper and Stanford, 2000; Benhamou et al., 2012; Mallet et al., 2012; Deister et al., 2013; Mastro et al., 2014; Abdi et al., 2015). The first round of network tuning used the original parameters of the STN model (Terman et al., 2002), which had autonomous firing rate around 2 spikes/s, but could not satisfy the constraints.
Since baseline firing rates depend on both cellular dynamics and synaptic conductances, we first adjusted the intrinsic dynamics of the cells, followed by synaptic parameter searches, iterating these two processes until our biophysical constraints were satisfied (see Fig. 1B). To further constrain our parameter sets, we analyzed the membrane potentials of the cells in the network and used only those that showed realistic nonpathologic voltage traces (see Fig. 1F), discarding, for instance, those which had any combinations of the following: very high peak voltages (>60 mV), hyperpolarized peak spike voltages, exaggerated AP half-widths, after-hyperpolarization potentials <−150 mV, RMPs depolarized >−50 mV, or more hyperpolarized than −85 mV. Final parameters are detailed in Table 2.
Cortex to STN stimulation
Our goal in this study was to examine both biophysical and architectural constraints that could give rise to the observed cortico-subthalamic theta power and spike rate modulations during response conflict. To this end, we first investigated cortical patterned activities that support increases in overall STN theta power, followed by how cortico-STN connectivity mediates theta power modulation in the context of conflict.
Cortico-STN drive
Cortical activity was modeled as spike trains activating STN AMPA and NMDA synapses. However, because of paucity of cortical spiking data in response conflict conditions, we examined two cortical spiking models (see Fig. 3A,B) that have been shown (Jones, 2016) to underlie frequency band power modulation.
The first cortical spiking profile, rhythmic single-spike drive (RSSD), consisted of single spike trains generated with fixed time periods, Tper, (rhythmic driving frequency,
The second spiking class consisted of a burst of spikes over a fixed duration EDur, referred to as single-burst event drive (SBED) (see Figs. 3B and 5B). The interspike intervals (ISIs) within the burst were drawn from a Gaussian distribution. In SBED, EDur is the same as the stimulation time Tstim. We vary spiking intensity, that is, the number of spikes in the SBED by controlling the means (IBIμ) and SDs (IBIσ) of the latter Gaussian distributions; thus, (
We extended the latter class to a rhythmic-burst events drive (RBED), where SBEDs are generated at fixed time intervals called interevent interval Δ (see Fig. 5B). In this case, the stimulation time,
Based on initial observations that replicated physiologically plausible synaptic responses and nonpathologic APs, we varied the NMDA conductances through two orders of magnitudes, with the values 5.208 × 10–7, 5.208 × 10–6, and 5.208 × 10–5 nS (specific somatic capacitance normalized to 1 μF/cm2), referred to as NMDALow, NMDAMid, and NMDAHigh respectively.
Cortico-STN architecture
Because we ultimately modeled response conflict in terms of simultaneous activation of multiple cortical populations representing different motor responses (see Simulating conflict), we also considered that these populations could in turn differentially project to STN subpopulations. Indeed, given the somatotopic organization of cortico-STN afferents (Nambu, 2011), we divided the STN population into four subpopulations and implemented two cortical feeds that probabilistically targeted two subpopulations. Since subpopulations conceptually represented somatotopic areas, cortical motor programs coding for mutually exclusive behaviors could thus preferentially innervate their own target STN subpopulations (see Fig. 5A). However, the degree to which such coding is completely segregated along cortical-STN “channels,” or whether crosstalk occurs to share cortical information, is unknown.
We thus parametrically varied the connection probability of a cortical feed to its own STN subpopulation. When a given cortical feed representing a single motor response connected exclusively to its own STN subpopulation, consequently segregating the cortico-STN communication, we call it a segregated topography (Fig. 5, showing only two subpopulations while four were simulated), with target probability ptar = 1 (all the STN units in the subpopulation are homogeneously driven by the same cortical feed) and nontarget probability pnon– tar = 0. A nonsegregated topography was defined by allowing crosstalk between cortical feeds and STN subpopulations by varying the value of ptar in the set (0.85, 0.55, and 0.25) under the constraint
ptar = 1 corresponds to a fully organized columnar topography, whereas ptar = 0.25 corresponds to a fully random one, as the cortical population was divided into four subpopulations. The probabilistic connectivities in the nonsegregated cases gave rise to STN units that received from both cortical feeds, termed conflict detectors (see Fig. 5A, right), since they could, at the single-cell level, listen to conflicting mutually exclusive motor programs. (Even in the segregated cases, however, it is still possible for the STN population as a whole to respond to conflict, given that multiple subpopulations would be coactive at once and could also indirectly interact via their reciprocal projections to GPe.)
Simulating conflict
To simulate conflict conditions, we activated both cortical feeds with temporal overlap (see Fig. 5B). To vary the levels of conflict, we varied the durations Edur of the feeds, as well as the delay δ and overlap ∩ between cortical feeds (see Fig. 5B).
Further, because many experiments report increases in cortical theta itself during response conflict (Cavanagh et al., 2014), which is Granger causal to STN theta (Zavala et al., 2014), we also tested the effects of rhythmic activity in cortical input, termed RBED (see Fig. 5B, right; see Cortical drives).
With these methods, we could thus investigate how STN biophysics, cortical-STN architecture, cortical dynamics, and their combinations impact STN theta.
Analysis
Custom Python scripts were written for all analyses and are available on request.
Spectral analysis
Spectral analysis was performed on simulated local field potentials (LFPs). Since we used only single-compartment units with no modeling of the extracellular medium and its resistivity, precluding any access to source and sink currents, and consequently also refrained from using arbitrary mappings to linearly convert from current to voltage to simulate LFP, and because LFP reflects mainly synaptic activities (Magill et al., 2004), the net synaptic currents served as a useful proxy for LFP. Thus, LFP was calculated as the mean squared of the sums of synaptic currents across all cells in the population of interest at every recorded time point, with units (nA)2. To contrast single STN unit LFP with no synaptic activity from GPeP, we alternatively used the net capacitative currents, for the panels in Figure 2 only. This was possible because the resonant frequencies agreed between both net capacitative and net synaptic currents in the lower bands (<80 Hz).
Spectral analysis was performed using the complex Morlet wavelet method, which was generated using Scipy and further normalized to ensure that the net input energy of the Morlet wavelet across various frequencies was controlled for. For every frequency in the range of 2-64 Hz (logarithmically distributed), we calculated the complex Morlet signal with 3.5 cycles and 217 (equivalent to ∼3.3-s-long kernel) points at a sampling period of 0.025 ms to generate a time frequency spectrogram.
To calculate theta energy, we integrated the theta power across stimulation time, which we report in units of (nA)2ms.
For efficiency, we opted for spectral domain multiplication using the Scipy fftconvolve routine after demeaning our time series data. The amplitudes of the resulting process were squared and treated as spectral power. To compute the power of a specific band, we average over all the frequencies that are members of that band at every time points.
Resonant frequency was defined as the frequency that showed highest power.
Spike statistics
All raster plots shown represent the exact recorded spike times (i.e., without binning); the latter were also used to compute the ISIs for every unit in the network.
To compute population spike counts, we first discretized the spiking data of each cell in fixed bin sizes of 5 ms, and aggregates in each bin were taken over all STN units that constituted a population of interest.
For cumulative spiking, we used the above spike counts for each subpopulation and cumulatively add them over succeeding time bins. To estimate the increase in spiking rate over the stimulation period, we divided the total number of spikes over that period by the total duration of the stimulation. For rhythmic bursts, the interburst intervals (IBIs) were also counted as part of the stimulation duration.
Statistical analysis
We performed simple two-tailed Welch's t test (SciPy Stats) between pairs of conditions, each with matched number of trials, for average theta power, average spiking rates and average theta-spiking, and we used a significance level of 0.05.
Results
Having tuned our model to reproduce empirical population spike statistics, we set out to explore the impact of synaptic conductances and cortical inputs on LFPs and spikes. We then proceeded to model “conflict” in cortical inputs, and how its effect on STN theta might be influenced by different architectural assumptions regarding cortico-STN topography.
The subthalamopallidal network requires cortically driven NMDA currents to show theta band resonance
We first confirmed that theta band activity did not emerge autonomously in the isolated STN unit or the connected network. Indeed, while the autonomous pacemaking activity of the isolated STN unit (i.e., ∼16 spikes/s when unconnected) eventually reproduced the empirical average firing rate of 11 spikes/s in conjunction with the inhibitory GPe layer (Fig. 1A,C–F), it did not show any theta band resonance. In this setting, the STN unit showed only higher frequency band activities (Fig. 2A), mostly reflecting harmonics of ∼16 Hz activity. To assess whether theta could emerge within the subthalamopallidal network (e.g., because of reciprocal GPe interactions; Fig. 1A), we averaged the net capacitative currents of all STN units and computed the spectrograms. Once again, no theta band power was observed (Fig. 2B). Because the network was probabilistically connected, we ran several (n = 15) instantiations; the absence of theta band activity was consistent across all such trials (data not shown; see also Fig. 1D).
A, STN-GPe network models with stochastic connectivities. Boxed numbers indicate connection probabilities. Ctx, Cortex (spike trains). B, Tuning process to achieve asynchronous baseline (no cortical drives) network spiking. C, Asynchronous baseline network spiking with STN raster and spike counts. Colors same as in A, time-aligned to cortical drive. D, Normalized spectral power of spike trains (in C). E, Top, ISI. Bottom, Spiking frequency counts of subpopulations in baseline network. F, Top, Autonomous cellular activities after tuning, in absence of any synaptic inputs. Bottom, Typical cellular activities in baseline network. G, Top, Normalized voltage-dependent synaptic conductances after tuning for AMPA, NMDA, and GABA_STN currents on STN, and GABA_GP on GPeP and GPeA. Bottom, Corresponding I–V characteristics, normalized to some peak currents.
We next tested whether and how cortical inputs to the STN could drive theta band rhythmicity within the network, as suggested empirically (Zavala et al., 2014). Because cortical rhythms are often expressed as transient bursts of activity (Jones, 2016), we considered multiple regimens of cortical drive: rhythmic single spike trains or as a burst of activity over a short duration (RSSD and SBED, respectively; see Cortical drives; for illustrations, see Fig. 3).
Because glutamatergic cortical drive could influence STN via both AMPA and NMDA currents, we isolated their contributions by blocking one or the other. These simulations revealed that AMPA currents are insufficient to generate theta band resonance even when cortical inputs are driven at 4 Hz (i.e., in the theta band; Fig. 2B). Notably, however, theta resonance emerged with the addition of NMDA synapses in this setting, and persisted with the removal of AMPA currents (Fig. 2B). Moreover, these results also held for SBED stimulations. As such, NMDA activity was necessary for theta band resonance in our model. We explore the biophysical mechanisms of this finding below, with particular focus on the NMDA currents, given their slower decay time courses (i.e., longer relaxation times) that could give rise to potent effects on theta.
A, Cellular activity of single unit after tuning. Top, Example is for a 4 Hz RSSD cortical drive. Other RSSD and SBED drives show similar dependence on NMDA for theta modulation. Left, Capacitative current (no synaptic components). Right, Spectrogram of capacitative current. B, Network activity by averaging capacitative currents across all STN units, under cortical inputs and in silico pharmacological blockade, aligned to cortical input onset.
Cortical bursts are filtered by NMDA dynamics to robustly elicit theta band activity
Having established that cortically driven NMDA currents promote STN theta band activity when stimulated in the theta range, we next examined whether theta was a resonant property of the driven network. A standard way to investigate resonance in nonlinear systems with feedback connections is to analyze its response to rhythmic (RSSD) external drives. We drove the network across a range of frequencies (2, 4, 8, 12.5, 20, and 32 Hz). We observed that, for mid and high levels of NMDA conductances, the STN showed resonance at the driving frequency (Fig. 3C), except at 32 Hz where these higher NMDA currents prevented repolarization and spiking. Interestingly, when the driving frequency was within the theta band (4-8 Hz), the resonant frequency matched the driving frequency for all levels of NMDA conductances. The dynamics of AMPA, NMDA, GABA, and net synaptic currents and spectrograms are shown for RSSD examples of 4 and 20 Hz in Figure 3A.
A, RSSD. B, SBED. Top, Cortical spike train. Middle, Synaptic currents (blue represents net; orange represents GABA; green represents NMDA; red represents AMPA) averaged across STN units in network. Bottom, Spectrogram. C, D, Resonant frequency (frequency with highest power) for (C) RSSD as in A and (D) SBED stimulations. C, D, Numbers indicate the corresponding cortical drives in A and B. Color intensity scales are different in each spectrogram.
We then tested the effects of cortical burst event patterns (SBED; see Cortical drives) as might occur during brief periods of response conflict (Isoda and Hikosaka, 2007, 2008) (without yet simulating conflict per se; i.e., these simulations were conducted with a single cortical feed). The durations and burst intensities of the events were varied to contrast a Low Stim protocol (10 ms duration and IBI of 9 ms, Fig. 3B, left) with a High Stim protocol (50 ms duration, 3 ms IBI, Fig. 3B, right). Despite this large variability in burst events, both simulations yielded theta band resonance. Confirming the low pass filtering effect of NMDA currents, the High Stim cases yielded lower resonant frequencies (Fig. 3D), because of larger NMDA peak currents' admitting longer relaxation times.
Notably, while both classes of cortical drives (RSSD and SBED) could produce theta band activities, these findings were differentially dependent on NMDA conductances (Fig. 3C,D). In the rhythmic driving mode, all levels of NMDA conductances facilitated theta resonance for RSSD in 4-8 Hz theta range, while for non-theta RSSD, higher NMDA conductances promoted resonance at the driving frequency. In contrast, burst event stimulations yielded theta band resonance over all three levels of NMDA and over a wide range of durations (albeit with lower resonant frequencies with longer event durations and higher burst spiking, as described above). In sum, cortical burst events robustly induced theta band resonance in the network, whereas rhythmic inputs required increasing NMDA conductance to do so.
Cortically evoked NMDA currents induce both STN burst spiking and silence periods (triphasic response) during theta resonance
What are the spiking characteristics of STN units during theta resonance? Consistent with empirical reports in both primates and rodents (Kitai and Deniau, 1981; Nambu et al., 2000; Hamani, 2004; Schmidt et al., 2013; Janssen et al., 2017; Pautrat et al., 2018) (Fig. 4E), STN units in our model could respond to cortical stimulation with either burst spiking (Fig. 4, left column) or with a triphasic spiking response characterized by initial synchronization, followed by a pause (defined as the average of the individual first ISIs that followed the first spiking evoked by cortical stimulation) and then finally bursts (Fig. 4, right column). The latter pattern was predominant for cortical spike trains that induced longer and larger NMDA currents, which were required for theta resonance. These simulations highlighted a critical role for this silent period in the triphasic response. Indeed, while one might expect that the silent period is induced by feedback inhibition from GPe, it occurred simultaneously with the time course of NMDA currents.
A-D, Left, Burst spiking response. Right, Triphasic (spike synchronization, silence, burst) spiking response, in STN subpopulation. A, STN raster and PSTH. Left ordinate, STN unit number. Right ordinate, spike count, bin size = 5 ms. B, Typical electrophysiological responses of stimulated units. C, Postsynaptic currents averaged across all STN units in subpopulation. D, Spectrograms (different color intensity scales). E, Published in vivo activities in rats from extracellular recordings: left, extracellular (Janssen et al., 2017); right, network (Schmidt et al., 2013) levels.
Given its role in theta resonance, we further examined the mechanism of this triphasic response. We analyzed individual STN membrane potential traces (Fig. 4B), and found that the pause in the triphasic response reflected a depolarization block induced by the rising amplitude of the NMDA current, followed by a membrane potential recovery period aligning with the decay time course of the NMDA current (Fig. 4C). Since the NMDA currents dominated the generation of the resonant theta frequency response, the duration of theta band resonance corresponded with the pause duration (Fig. 4D). Only after NMDA currents decayed sufficiently, the GPe-mediated GABAergic currents became effective, leading to subthreshold oscillations in the membrane potential, followed by burst firing, completing the triphasic response.
To further understand the potentially separable roles of NMDA and GABA currents in this cortically driven triphasic response, we varied both of their conductances. We first confirmed that increases in GPe GABAergic conductance increased the pause durations (Hamani, 2004) and reduced the burst activities, without affecting initial spiking (data not shown). Moreover, larger NMDA conductances also prolonged the silence period. Given their opposing effects on postsynaptic potentials, that both contributed to lengthen the pause in spiking was rather counterintuitive; nevertheless, different mechanisms underlie these similar effects.
Finally, the post-silence bursts are thought to reflect rebound bursting because of T-type calcium currents activated by removal of GPe-mediated hyperpolarization (Hamani, 2004; Bevan et al., 2006, 2007; Magill et al., 2006). Our investigations contribute an additional mechanistic explanation for such triphasic spiking activity. Congruent with the extant interpretations, it is initiated by excitatory AMPA currents that promote initial spiking; the onset of pause is triggered by synchronous recruitment of GABA activities from GPe. However, our model suggests that the pause is sustained beyond this inhibition because of the depolarization block induced by NMDA currents that weaken the effects of K+ currents, thereby preventing the hyperpolarization necessary for AP generation. While this prediction needs to be experimentally tested, it accords with previous studies (Wilson et al., 2004; Barraza et al., 2009) that report depolarization block in the STN and their link to potassium currents. We observed that as long as cortical spiking persisted, the membrane potential was held at ∼−20 mV while the NMDA current accumulated. Once cortical spiking ceased, the NMDA current started to relax with a slow time course, prolonging the pause. Since the latter dynamics did not feature hyperpolarization, rebound bursting was because of depolarization activation of L-type calcium currents. See Beurrier et al. (2001) for alternative mechanisms underlying the silence period following high-frequency stimulations, and Magariños-Ascone et al. (2002), Perlmutter and Mink (2006), Chiken and Nambu (2016), and Hamani et al. (2017), for related “depolarization block” mechanisms involved in STN deep brain stimulation in PD.
Theta power is strongly modulated by conflict and cortical to STN topography
With a clearer understanding of how NMDA contributes to theta power modulation, we tested how such modulation may be additionally influenced by cortico-STN topography (see Cortico-STN architecture; Fig. 5A), and response conflict (Fig. 5B). To simulate response conflict, we included two populations of presumed cortical motor units (representing mutually incompatible responses), based on a definition of conflict as Hopfield energy (Botvinick et al., 2001; Yeung et al., 2004; Frank, 2006). In the low-conflict condition, only one of these populations was active, but in a high-conflict condition, they were coactive with increasing temporal overlap (Fig. 5B). We varied topography of the cortical inputs, ranging from completely segregated (each cortical population projects singly to its target STN subpopulation) to fully random connections (each STN unit is equally likely to receive inputs from any cortical population; Fig. 5A). The number of STN units that receive inputs from more than one population increases as the topography becomes more random and less segregated. We posited that such STN units could act as “conflict detectors” because they would detect the coactivation of multiple cortical populations and allow us to assess whether they would impact theta power in accordance with the literature (see Introduction). Cortical activity was simulated as burst drives in the two populations as would occur during action selection (Fig. 5B), either as single events (SBED) or repeated “rhythmic” burst events (RBED), the latter representing vacillating conflict within the cortex, consistent with the literature on frontal theta and conflict.
A, Cortico-subthalamic network architecture. Left, Segregated information flow. Cortical subpopulations connect to specific STN subpopulations, precluding direct information sharing, showing only two subpopulations, but four were used in the simulations. Middle, Nonsegregated information flow. Cortical subpopulations connect primarily to their corresponding STN subpopulations with high probability while allowing connections to other STN subpopulations with lower probability. Right, Different STN classes based on stochastic cortical-STN connectivities. Classes: MainStim, STN units that receive only from their own cortical subpopulation; OtherStims, STN units that receive from only one cortical subpopulation that is not their own; NoStims, STN units that do not receive any cortical inputs; AllStims, STN units that receive from more than one cortical inputs, acting as conflict detectors. B, Simulating conflict condition. Top, No conflict; only one cortical subpopulation active. Bottom, Conflict; coactivation of two cortical subpopulations, with offset δ. Left, Stimulation schematics. Middle, SBED. Right, RBED.
Notably, overall theta power in the network increased with lower segregation, that is, with higher number of conflict detectors (Fig. 6A). An example of such high conflict with a delay between two cortical events (δ) of 5 ms is shown in Figure 6B (two-sided Welch's t test, p < 0.05). This result is not simply because of more input; indeed, we maintained the total number of inputs to the STN to be equal across all levels of segregation. Moreover, in the “no conflict” condition, theta was actually maximized by the highest segregation level, thus suggesting that topography favoring conflict detectors enhances theta only when conflict is actually present.
A, Average LFP theta power across all STN units that received cortical inputs. Left, No conflict. Right, Conflict. B, Average LFP theta power across different STN classes. Subpopulations first and second received from either the first or the second cortical subpopulations only. Both, STN conflict detector units received from both cortical drives; 1st + 2nd, the (linear) average power from subpopulations first and second lumped together. C, Average synaptic currents. Left, AMPA. Right, NMDA, across STN detector units. Color coding same as in B. D, Effects of cortical stimulation offset (δ) on average theta. Left, Power. Right, Energy. E, Contribution of glutamatergic currents. Left, NMDA. Right, AMPA, toward theta power. Top, δ = 5 ms. Bottom, δ = 10 ms. Broken lines indicate time points 0, 5, and 10 ms. F, Nonlinear effects of cortical stimulation overlaps and durations, and segregation on the relative charge contributions of NMDA and AMPA currents observed during theta activity.
NMDA currents enhance conflict-induced theta via supralinear summation
These results recapitulate the dominant finding in the human STN LFP literature (Cavanagh et al., 2011, 2014; Zavala et al., 2013, 2014, 2016; Herz et al., 2016), whereby response conflict increases STN theta power. Our network model allows us to investigate whether such theta modulation was preferentially increased in conflict detectors. We thus computed the average theta power in separate STN subpopulations: those receiving inputs from either of the individual cortical drives (first or second) and those receiving from both (1st + 2nd). We found that theta power was indeed preferentially elevated in conflict detectors at all segregation levels. This finding stems from a supralinear summation afforded by the biophysical properties of the STN network: theta power in conflict detectors was greater than the sum of power in two populations that each received from only one cortical feed (see, e.g., Fig. 6B).
To understand the underlying electrophysiological mechanisms for such supralinear summation, we analyzed the different glutamatergic currents, AMPA and NMDA (Fig. 6C). We found that NMDA currents generated in the conflict detectors were indeed greater than the sum of the currents from the individual subpopulations receiving single cortical inputs. This pattern was not seen in AMPA currents (indeed, it was somewhat reversed, with smaller currents in the conflict detectors compared with segregated populations). This analysis suggests that response conflict increases STN theta via supralinearity mediated by NMDA currents in conflict detectors.
To further investigate how the NMDA-based dynamics of the conflict detectors interact with cortical spiking, we manipulated the durations (10-125 ms) and temporal offsets (delays) δ, (10-37 ms) (Fig. 5B) of the cortical stimulations, and consequently, the temporal overlap (i.e., conflict as defined above). We found that theta power was maximized for the lowest delays (highest overlap) in the unsegregated network (most conflict detectors). In contrast, a segregation level of 55% was more apt at theta maximization for higher delays (δ ≥ 10 ms). Figure 6D shows theta power and theta energy for SBED durations of 10 ms; these results are robust for SBED durations of up to 80 ms.
To understand how delay, δ, is critical in encoding conflict and impacting theta, and how it interacts with segregation, we analyzed the glutamatergic currents, and found that again, NMDA currents but not AMPA currents, are directly modulated by δ (Fig. 6E). In particular, the lower the delay, the higher the opportunity for NMDA currents from the different SBEDs to integrate supralinearly, thus yielding higher theta with lowest segregation. AMPA currents show no modulation with respect to either delays or segregation levels. This supports the idea that high-conflict condition is encoded by the relative timing of incoming cortical signals, and mechanistically realized through the nonlinear NMDA dynamics.
It was surprising that only NMDA and not AMPA would show modulations with respect to segregation, despite being driven by identical cortical signals. Since this is because of their respective excitatory postsynaptic current rise and decay kinetics and the membrane potential synaptic conductance modulations (see Synaptic tuning), we characterized their differential effects by computing the ratio of overall charge transfer of NMDA to AMPA, qNMDA:qAMPA (Fig. 6F). At the lowest delay of 5 ms, the ratio increased with cortical signal overlap, ∩ (see Fig. 5B). Moreover, when overlap was low, there was a positive correlation between the ratio and number of conflict detectors (hence, inverse relationship with segregation levels). However, at higher delays, δ ≥ 10 ms, mid-level segregation was better at reflecting the contrasting contributions. Since combinations of δ and SBED durations would lead to conditions of zero overlap ∩ = 0 ms, we repeated the same analyses against the total duration time of stimulations that were provided to the units. The results showed similar trends as above, further confirming that mid-level segregation provides a wider dynamical range as a signal integrator when the delays are large.
These results are notable in that the literature implicates higher theta powers being linked to higher response time levels during response conflict, while this pattern is not seen, and indeed is often reversed, in low-conflict situations (Cavanagh et al., 2011; Herz et al., 2016). We return to this issue with a plausible explanation in the Discussion.
Rhythmic burst events at theta frequency maximize STN spiking and theta power
Above, we observed that STN theta power could be expressed either when cortical inputs themselves oscillate in the theta range (RSSD), or when they consist of single burst events (SBED), especially with conflict. However, we also noted that higher theta power was associated with longer pauses in the triphasic spiking response. This latter result is counterintuitive, as several studies have reported increased STN spike rates during conflict, which correlate with larger decision thresholds (Isoda and Hikosaka, 2008; Zaghloul et al., 2012; Herz et al., 2016). To investigate a possible reconciliation of these observations, we first noted that, in our model, STN bursting occurred after NMDA currents decayed, which also signaled the end of theta band resonance. We then assessed the impact of prolonged response conflict in which rhythmic and burst event modes are combined (RBED), representing a situation in which conflict persists for multiple cycles in cortex (Fig. 5B), consistent with the cortical literature (Bartoli et al., 2018; Zeng et al., 2021). We tested this notion by presenting repetitions of high stimulation burst events (50 ms duration, 3 ms IBI, 12 burst events), where the events are repeated with different intervals Δ tested at 75, 125, 200, 350, and 450 ms. Notably, time periods from 125 through 250 ms fall in the theta range (4-8 Hz) period.
We first tested this rhythmic event (RBED) spiking protocol in the “no conflict” condition as above and observed that, unlike in the previous case, theta power did not differ significantly across segregation levels (results not shown). Next, we presented two simultaneous rhythmic event drives to examine the influence of conflict (Fig. 7A, top). Theta power was higher for every IBI compared with the “no conflict” condition, as expected (data not shown but were similar to the contrast in Fig. 6A). The instantaneous theta power was characterized by a prototypical initial peak, followed by an oscillatory period, the amplitude of which increased with IEI (Fig. 7A). As previously, theta amplitudes were strongest for the lowest segregation levels, across all IEI, Δ. During the oscillatory instantaneous theta period, we observed sustained theta amplitude levels that peaked when Δ was in theta range time period (125, 250 ms). When Δ was higher than that, instantaneous theta would rise but then decay to 0, leading to lower overall sustained power.
Effects of interevent intervals (Δ) in RBED cortical stimulations and segregation on STN conflict detector units. A, Top, Average LFP theta power. Bottom, Cumulative average spike counts, during and after cortically induced conflict. Overlayed tables represent the net spiking rates during (=) and after (>) stimulation, and ratios of spiking rates during/before (=/<) and during/after (=/>) stimulation periods. B, Box plots represent distributions (orange line indicates median) and trial averages (green triangles represent mean) of conflict detectors during RBED stimulation period. Top, Average theta energy. Middle, Average spiking. Bottom, Average theta spiking (theta × spiking). Within each segregation level, effects of Δ on theta energy, spiking, and theta-spiking are significant (p < 0.05) between theta-IEI 125, 250 and non-theta-IEI groups 75, 350, 450. C, Average cortically driven postsynaptic currents of conflict detectors, shown here for the lowest segregation level 25%. Colors represent IEI. NMDA, but not AMPA, responses show aliasing because of nonlinear interactions.
To better characterize these nonlinear dynamical effects, we computed the overall theta energy (trial-averaged instantaneous theta power integrated over the stimulation period). First, as noted above, theta energy was highest in the lowest segregation level across all Δ and generally decreased with higher segregation (Fig. 7B, top), although the effect of segregation was not significant (consistent with the results in Fig. 6F where we observed that longer stimulation durations would eventually blur the segregation distinctions). Second, within each segregation level, Δ in theta range (125-250 ms) produced largest overall theta energy compared with non-theta ranges (75, 350, 450 ms). The results were significant (two-tailed Welch's t test, p < 0.05) between these two range groups. These were because the theta-range promoted more elevated sustained theta, while non-theta ranges produced lower sustained theta (Δ = 75 ms) or instantaneous theta dropping to zero (Δ = 350, 450 ms) that decreased their contribution to theta rhythmicity.
Finally, since STN spiking increases during response inhibition (Isoda and Hikosaka, 2008; Schmidt et al., 2013; Fife et al., 2017) and with conflict-induced elevations in decision thresholds (Herz et al., 2016), we tested the effects of IEI on spiking. We found that average spike rates of conflict detectors during cortical stimulation were maximized, and ratios of spiking rates during conflict with respect to pre- (≤) and post- (≥) conflict conditions were highest, for Δ lying in the theta period range (125-250 ms) across all segregation levels (Fig. 7A, bottom, B, middle), with results being significant between the two range groups (two-tailed Welch's t test, p < 0.05). While there was a clear monotonic trend for theta energy as regards segregation, we observed more variable effects on spiking among segregation levels. Indeed, this is consistent with the above observations that spiking and theta power are not always positively correlated, the more so as NMDA conductance increases.
We thus defined a new metric, theta-spiking, by taking the product of theta energy and spiking rates (Fig. 7B, bottom), to understand under what conditions both theta and spiking rates would be concurrently maximized, which would be consistent with response conflict neural dynamics. First, we observed that theta-spiking was maximized when Δ was in theta range (125-250 ms), with significant differences (two-tailed Welch's t test, p < 0.05) between the two range groups as observed for both theta energy and spiking rates separately. Moreover, segregation showed a clear monotonic trend, with lower segregation level of 25% maximizing theta-spiking while 85% minimizing it.
To understand the underlying biophysical mechanisms mediating these processes, and since theta power here is a direct monotonic reflection of the underlying cortically driven postsynaptic currents, we again focus on the net glutamatergic responses (Fig. 7C). Notably, for NMDA kinetics, there was a nonmonotonic relationship with Δ, whereby its dynamics favor IEI lying in theta range (125-250 ms). IEIs lower than that lead to aliasing of cortical NMDA currents, whereas those larger than the theta range exhibit no aliasing but have extended periods of silence. Δ = 200 ms lies in an optimal range that not only prevents aliasing, which enhances theta power, but also allows the NMDA currents to decay enough to basal levels to eventually promote higher spiking rates. In contrast, AMPAergic response to burst spiking remained invariant with IEI, both in trajectory and durations (∼80 ms), owing to their fast rise and decay kinetics in conjunction with the voltage-dependent synaptic conductance modulation (Fig. 1G; see Synaptic tuning).
In sum, we observe that high conflict, high overlap and the presence of conflict detectors maximize theta. NMDA dynamics provide a filtering effect of the stochastic network spiking, thus robustly generating theta in these cases. In contrast, STN spiking is more susceptible to ongoing network stochastics, occurring mainly when NMDA currents are relatively low at the cellular level, but nevertheless maximized when cortical inputs are in the theta range. These results are consistent with observations that cortical theta is Granger causal to STN theta (Zavala et al., 2014), while also accounting for STN-GPe cellular and network dynamics that mediate changes in spiking.
Discussion
We provide a novel theoretical framework for mechanistically interpreting the growing body of evidence linking STN theta oscillations to conflict and decision threshold adjustment (Cavanagh et al., 2011; Green et al., 2013; Zavala et al., 2013, 2014, 2016, 2017; Frank et al., 2015; Herz et al., 2016, 2017; Ghahremani et al., 2018; Kelley et al., 2018; Wessel et al., 2019;). Such findings accord with existing neural network models of BG, in which increased STN spike rates induce a transient increase in decision threshold during response conflict, by modulating thalamocortical activity (Frank, 2006; Ratcliff and Frank 2012; Wiecki and Frank, 2013; see also Bogacz and Gurney, 2007). Our biophysical explorations of STN-GPe circuitry lend insight into the underlying mechanisms, leading to the following novel and testable predictions.
First, our model suggests that STN theta power does not spontaneously emerge but rather requires cortical input and is strongly dependent on NMDA currents. The dynamics of cortical inputs were also critical: STN units exhibited theta band resonance if the cortical inputs themselves oscillate in theta (Zavala et al., 2013, 2014; Cavanagh et al., 2014; Wessel et al., 2019), but theta also robustly emerged in response to cortical burst events at a broad range of intensities and durations, particularly when such events overlapped (mimicking conflict). Moreover, both theta power and STN spiking were elevated and prolonged during rhythmic burst events, representing a state of conflict in which cortical motor plans vacillate in the theta range. These results accord with findings that individual STN neurons exhibit increased spike rates to conflicting evidence (Isoda and Hikosaka, 2008; Zaghloul et al., 2012), and that cortical theta is Granger causal to STN theta (Zavala et al., 2016). Finally, theta band resonance concomitant with spiking was also strongly modulated by architectural constraints, with maximal response when cortical inputs provided divergent connectivity to multiple STN subpopulations. Analysis of the underlying mechanisms of such effects revealed an NMDA-dependent supralinear response in STN “conflict detector” units.
NMDA mechanisms
The slow time course of NMDA currents, biophysically constrained from the literature, was responsible for the modulation of theta. Theta band resonance emerged even in response to single burst events, and magnified during rhythmic bursts, because of the time period of the envelope of the decaying phase of the NMDA current. We did not tune any parameters to reproduce these experimental observations; they emerged from the biophysical constraints. Notably, while theta band resonance is often described in terms of low frequency “oscillations,” our investigations align with biophysical models demonstrating how resonance can also emerge from burst events and not just rhythmic activities (Jones, 2016).
On the applied end, such findings imply that it may be possible to enhance STN function in tasks that require cognitive control by targeting local NMDA receptors, particularly the NR2D subtype, which are most prevalent in STN and GPe (Callahan et al., 2020; Yi et al., 2020). Indeed, these authors have recently developed a positive allosteric modulator of NR2D, and showed that it effectively reduced premature responding in a conflict task in rats, precisely the phenotype that is elevated with STN lesions (Baunez and Robbins, 1997; Baunez et al., 2001), and which motivated the theory that STN implements a “hold your horses” mechanism (Frank, 2006). Our biophysical model provides a rationale for why NMDA currents would influence such function, and also suggests that yet stronger doses could reverse the effect (because of depolarization block). Thus, in principle, different doses of such an agent could be useful for both disorders of impulsivity (too low decision threshold) and compulsivity (too high).
Cortico-STN communication
We found that STN theta was particularly elevated when multiple cortical populations targeted individual STN units and overlapped in time (i.e., conflict). The primate STN is topographically organized by functional domains with broadly defined boundaries, but there also exist convergence zones (Nambu, 2011; Haynes and Haber, 2013; Kita et al., 2014; Alkemade et al., 2015). While our findings implicate STN units as conflict detectors, STN also receives direct projections from dorsomedial frontal areas (pre-SMA and anterior cingulate), which themselves resonate in theta during conflict (Aron et al., 2007; Zavala et al., 2014; Frank et al., 2015). Future work should explore separable roles of these cortical inputs, but our rhythmic burst simulations suggest that such cortical theta inputs may amplify STN theta and spiking. These interactions may also be related to observations that high-frequency stimulation in PD (through STN excitation or silencing) can be effective by disrupting low-frequency signals but preserving high-frequency ones (Garcia et al., 2005).
Reconciling the theta-response time conundrum
While many studies have replicated the finding that low-frequency oscillations are related to increases in response time and decision threshold, it is noteworthy that all of these findings were specific to high-conflict task conditions (Cavanagh et al., 2011; Green et al., 2013; Zavala et al., 2013, 2014, 2016, 2017; Frank et al., 2015; Herz et al., 2016, 2017; Kelley et al., 2018; Wessel et al., 2019). If frontal and STN theta band powers were simply a “read out” of experienced conflict, then one should be able to predict response time from theta power, regardless of task condition. However, empirical findings have shown the opposite: within low-conflict conditions, larger frontal and STN theta powers are related to reduced response time (Cavanagh et al., 2011; Herz et al., 2016). Nevertheless, evidence abounds across species that STN is needed for stopping action rather than facilitating it (Baunez et al., 2001; Aron, 2007; Aron et al., 2007; Isoda and Hikosaka, 2008; Jahfari et al., 2012; Fife et al., 2017; Jahfari et al., 2019; Wessel et al., 2019). Our model provides plausible mechanistic reconciliations of this observation.
First, we showed that theta was observed primarily in conflict detector units. In low-conflict situations, there was still an overall increase in theta with increased cortical intensity (driven by NMDA currents in individual STN populations), but that increased theta at the population level does not necessarily involve increased spiking across the STN population. Indeed, empirical data show that STN spike rates increase in conflict conditions (Isoda and Hikosaka, 2008; Zaghloul et al., 2012); moreover, disrupting STN function only impairs response inhibition during high-conflict or surprising situations (Frank et al., 2007; Cavanagh et al., 2011; Fife et al., 2017; Ghahremani et al., 2018). Notably, in our model, STN theta is related to spike output in opposite directions, depending on the period of time after stimulus onset because of the triphasic response. Early during a decision, theta rises during the silent period, and during low-conflict trials, the large evidence in cortical input may be sufficient to induce fast responding. In high-conflict conditions, choices will not yet have reached threshold; the following burst spiking period, which is amplified with rhythmic cortical bursts, induces an increase in effective decision threshold. This interpretation is consistent with neural network models in which STN increases and then collapses, leading to dynamic decision thresholds over the course of choice (Ratcliff and Frank, 2012). Moreover, such a dynamic threshold is normative specifically in tasks that involve a mixture of low and high-conflict trials (Malhotra et al., 2018). Thus, STN theta power may speed up choice during low-conflict situations simply because spiking is actually reduced early during a choice process. Finally, it is also possible that BG is not involved as much in action selection when the correct action is sufficiently salient (Brown and Marsden, 1998; François-Brosseau et al., 2009; Cockburn et al., 2014).
In conclusion, in these studies, we have provided both biophysical and architectural mechanisms associated with theta band power modulation by focusing on cortically evoked STN-GPe interactions. We have deliberately not included other parts of the BG network, including the striatum, to focus on the mechanisms of cortically driven theta power in STN which appear to be critical for modulating response time during conflict tasks. This does not imply that other parts of the network are irrelevant, and indeed it would be interesting to investigate the roles of the D1 and D2 striatal pathways on GPe and STN communication, given their strong implication in motivated behavior through cortical innervations that might also converge on the STN. For a more comprehensive picture, the model could be extended to include striatal effects to investigate how the triphasic spiking responses are differentially modulated by the D1 versus D2 tug-of-war, with particular emphasis on the arkypallidal GPe subpopulation.
A paucity of studies characterizing cortical spiking during response conflict has also limited our explorations to just a few potential cortical spiking profiles. We would expect the BG-thalamocortical loop to have more dynamical effects on cortico-STN communication, which we did not account for. Moreover, although synaptic activity would monotonically correlate with LFP, we could not consider spatial filtering effects because we neither simulated the oblong arrangement in the STN nor its detailed cellular geometries, which might affect the resonant frequencies.
Footnotes
This work was supported by the Fulbright scholarship, the Open Graduate Education program at Brown University, the Brown Institute for Brain Science Research Graduate Award (supported by the Dr. Daniel C. Cooper Graduate Award, the Charles A. Dana Graduate Fellowship Fund, and the Department of Psychiatry and Human Behavior), and National Institute of Mental Health Grant R01MH115905-01, National Science Foundation Proposal 1460604, and National Institutes of Health Grants RO1 MH106174 and P20GM103645. We thank Maxwell Sherman and Shane Lee for coding support; Prof. Jonathan E. Rubin, Prof. Nancy J. Kopell, Prof. Barry W. Connors, and Assoc. Prof. Wael F. Asaad for comments; and the Brown University Center for Computation and Vision, and the Neuroscience Gateway for computing resources.
The authors declare no competing financial interests.
- Correspondence should be addressed to Michael J. Frank at Michael_Frank{at}brown.edu
