QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

HOME  |   SEARCH  |   ARCHIVE  |   SUBSCRIBE  |   CONTACT  |   HELP

The Journal of Neuroscience, December 13, 2006, 26(50):12921-12942; doi:10.1523/JNEUROSCI.3486-06.2006

This Article Abstract Full Text (PDF) Submit an eLetter Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by Humphries, M. D. Articles by Gurney, K. N. Search for Related Content PubMed PubMed Citation Articles by Humphries, M. D. Articles by Gurney, K. N.

 Previous Article  |  Next Article 

Behavioral/Systems/Cognitive
A Physiologically Plausible Model of Action Selection and Oscillatory Activity in the Basal Ganglia

Mark D. Humphries, Robert D. Stewart, and Kevin N. Gurney

Adaptive Behaviour Research Group, Department of Psychology, University of Sheffield, Sheffield, S10 2TP, United Kingdom

	Abstract

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The basal ganglia (BG) have long been implicated in both motor function and dysfunction. It has been proposed that the BG form a centralized action selection circuit, resolving conflict between multiple neural systems competing for access to the final common motor pathway. We present a new spiking neuron model of the BG circuitry to test this proposal, incorporating all major features and many physiologically plausible details. We include the following: effects of dopamine in the subthalamic nucleus (STN) and globus pallidus (GP), transmission delays between neurons, and specific distributions of synaptic inputs over dendrites. All main parameters were derived from experimental studies. We find that the BG circuitry supports motor program selection and switching, which deteriorates under dopamine-depleted and dopamine-excessive conditions in a manner consistent with some pathologies associated with those dopamine states. We also validated the model against data describing oscillatory properties of BG. We find that the same model displayed detailed features of both -band (30–80 Hz) and slow (1 Hz) oscillatory phenomena reported by Brown et al. (2002) and Magill et al. (2001), respectively. Only the parameters required to mimic experimental conditions (e.g., anesthetic) or manipulations (e.g., lesions) were changed. From the results, we derive the following novel predictions about the STN–GP feedback loop: (1) the loop is functionally decoupled by tonic dopamine under normal conditions and recoupled by dopamine depletion; (2) the loop does not show pacemaking activity under normal conditions in vivo (but does after combined dopamine depletion and cortical lesion); (3) the loop has a resonant frequency in the -band.

Key words: neural network; action selection; dopamine; oscillations; subthalamic nucleus; globus pallidus

	Introduction

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Work from numerous fields is converging on the hypothesis that one function of the basal ganglia (BG) is the selection of motor programs (Chevalier and Deniau, 1990; Mink and Thach, 1993; Hikosaka et al., 2000). A more specific hypothesis is that the BG are a critical neural substrate in the vertebrate action selection system, resolving conflict between multiple neural command centers trying to access the final common motor pathway (Prescott et al., 1999; Redgrave et al., 1999). A range of computational modeling work, from our laboratory (Gurney et al., 2001b; Humphries and Gurney, 2002) and others (Beiser and Houk, 1998; Frank, 2005), has demonstrated results consistent with the selection hypotheses. Previous models have abstracted much biological detail by using either population-level representations or mean firing rate elements throughout. As such, no existing model makes contact with the substantial body of work detailing the firing properties of BG neurons (Bevan et al., 2002; Boraud et al., 2005). To address this, we recast the Gurney et al. (2001b) model with populations of more realistic "spiking neurons."

This strategy allowed us to study the prevalent issue of the oscillations that characterize neural activity of the BG in both motor disorders and motor control (Boraud et al., 2005). Many oscillatory phenomena occur in the subthalamic nucleus (STN) and the globus pallidus (GP), whose neurons are reciprocally connected (Bevan et al., 2002). Prominent -band (30–80 Hz) oscillations have been found in local-field potential (LFP) recordings in the BG of both awake and anesthetized animals, including humans (Brown et al., 2001; Berke et al., 2004; Magill et al., 2004b). Their absence in the BG of Parkinson's disease patients (Brown et al., 2001) and increase in power during motor activity (Brown et al., 2002; Cassidy et al., 2002) suggests that -band oscillations are intimately related to normal motor behavior (MacKay, 1997). However, despite the anatomically identified feedback loop formed between STN and GP, the two nuclei seem decoupled in the healthy BG (Magill et al., 2000; Raz et al., 2000; Urbain et al., 2000; Baufreton et al., 2005). For example, slow (<1 Hz) cortical oscillatory activity is reflected in STN neuron output but not transmitted to GP neurons, except when dopamine has been depleted, yet, after removing both cortex and dopamine, some residual slow oscillatory activity remains in both nuclei (Magill et al., 2001). Such data seem at odds with hypotheses for the role of the STN–GP loop in motor program selection (Berns and Sejnowski, 1996; Bevan et al., 2002; Rubchinsky et al., 2003).

The purpose of the present study was to generate a new spiking neuron model of the principal intrinsic connections of the motor BG, test its capacity as an action selection mechanism, study its oscillatory properties, and assess them against relevant biological datasets. We therefore validated the model against experimental studies by Brown et al. (2002) and Magill et al. (2001) on -band and slow oscillatory phenomena, respectively. When possible, parameter values were derived from experimental studies. Furthermore, in all simulations, the majority of the parameter values were unchanged: only those required to emulate the experimental conditions were changed. We also sought to use the results of the model to propose underlying neural mechanisms for the oscillatory phenomena.

Parts of this work have been published previously in abstract form (Humphries et al., 2004; Gurney et al., 2005).

	Materials and Methods

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The model
We base our model on the biologically informed functional anatomy of the BG proposed by Gurney et al. (2001a), which is depicted in Figure 1. When possible, all anatomical and physiological data used to constrain the model have been derived from rat-based studies, and so we use rat-brain terminology throughout. The nuclei and connectivity we incorporate follow current views of BG organization (for example, see Bolam et al., 2000), which we itemize as follows.

View larger version (30K): [in this window] [in a new window]   Figure 1. Functional anatomy of the BG. Excitatory glutamatergic cortical input drives cells in the two input nuclei of the BG, the striatum and the STN. The GABAergic striatal projection neurons can be subdivided into two populations on the basis of the dominant dopamine receptor type (D1 or D2 type) that modulates their firing (see Materials and Methods). The D1-dominant cells predominantly send inhibitory projections to the output nuclei, the SNr (shown here) and the EP nucleus (the model presented includes only the SNr, because the two nuclei are similar with respect to their intra-BG connections). These GABAergic cells in turn send tonic inhibitory output to their targets in the thalamus and brainstem. It is the removal of this tonic inhibition that is thought to signal the selection of an action; hence, the combination of striatal D1 cells and the output nuclei are termed the selection pathway. The inset box shows the proposed off-center, on-surround network formed by this pathway, with hypothetical activity of the channel populations. The D2-dominant cells of the striatum send inhibitory projections to the GABAergic cells of the GP; the GP cells in turn send inhibitory projections to the STN and SNr/EP and are a source of control signals for the basal ganglia (hence, this circuit is termed the control pathway). Glutamatergic cells of the STN send excitatory projections to the GP and the SNr/EP, contributing to both the tonic activity of the SNr/EP cells and the activity within the control pathway.

 
(1) The striatum is considered the principal input nucleus of the BG, receiving extensive input from most cortical regions (Gerfen and Wilson, 1996; Glynn and Ahmad, 2002) and thalamic nuclei (Smith et al., 2004). GABAergic projection neurons comprise 95% of the neuron population in striatum (Gerfen and Wilson, 1996).

(2) The STN forms the other primary input nucleus of the BG and also receives input from cortex and some thalamic nuclei. Only one neuron type is known to exist in STN, and all cells are glutamatergic and project outside of the nucleus (Temel et al., 2005).

(3) Both striatum and STN projection neurons send axons to GP and the output nuclei of the BG, the substantia nigra pars reticulata (SNr), and the entopeduncular (EP) nucleus (Bolam et al., 2000).

(4) The GP neurons are GABAergic and project to the STN (so that STN and GP form a closed loop) and to the output nuclei (Smith et al., 1998).

(5) The GABAergic cells of the output nuclei contact numerous thalamic nuclei and brainstem structures (Bolam et al., 2000).

(6) The partition of the striatum into two projection neuron populations based on their dominant dopamine receptor type (D₁ or D₂ type) (Bolam et al., 2000) is maintained in the current model. Significant evidence exists for the colocalization of these receptors in some or all of the projection neurons (Surmeier et al., 1996; Aizman et al., 2000). However, many converging lines of evidence from electrophysiology (Gonon, 1997; Kaneda et al., 2002; West and Grace, 2002), microdialysis (O'Connor, 1998), mRNA transcription (Gerfen et al., 1990; Murer et al., 2000), and lesion (Sano et al., 2003) studies suggest a functional split between D₁- and D₂-dominant projection neurons and, furthermore, that the D₁-dominant neurons project to SNr/EP and the D₂-dominant neurons project to GP.

(7) A central concept governing the organization of BG connectivity is the existence of parallel anatomical loops running throughout the BG nuclei (Alexander et al., 1986; Middleton and Strick, 2000). Furthermore, a distinction can be made between macroscopic and microscopic channels. Macroscopic channels are formed by closed send-and-return loops running in parallel, each originating from an identified cortical area, passing through the BG, and returning to the originating cortical area via thalamus (Alexander et al., 1986). [These loops are also open in the sense that projections from different cortical areas converge on the same locations in striatum (Romanelli et al., 2005) and that there may be a hierarchy of the loops created by interconnections within the BG (Joel and Weiner, 2000; Haber, 2003)]. Recent trans-neuronal tracing data are consistent with at least 10 such channels within the BG output nuclei (Middleton and Strick, 2000), although it is not clear to what extent they are sustained within the BG intrinsic circuitry. Some authors are content to divide the BG intrinsic circuitry into three parallel domains (Joel and Weiner, 2000), following a rough subdivision of the striatum into motor, associative, and limbic territories, based on the send-and-return loops of whole cortical regions (respectively, motor cortices and the division of prefrontal cortices into associative and limbic areas).

Microscopic channels are discrete parallel loops running within a macroscopic channel. For example, the somatotopic map found within the striatal motor territory is maintained throughout the BG intrinsic circuitry, such that there are separate channels for arm, leg, and face representations (Alexander and Crutcher, 1990; Romanelli et al., 2005). Similar topographic maps have been proposed for the other macroscopic channels (Alexander and Crutcher, 1990). Moreover, within these limb representations, there may be additional discrete channels corresponding to particular movements, demonstrated in striatum by microstimulation (Alexander and DeLong, 1985) and markers for metabolic activity during behavior (Brown and Sharp, 1995), and/or to detailed body map locations (Brown et al., 1998). We model here such microscopic channels, each channel representing a different putative action (Redgrave et al., 1999; Gurney et al., 2001a). Anatomically, the channels in striatum and STN are defined by the regional extent of the topographically fragmented input from distinct representations in cortex to these areas and the channels in GP and the output nuclei by their corresponding striatal afferents.

(8) We proposed that the amount of activation in each channel at the level of the striatum effectively represents the relative "salience" or urgency of the action currently represented at that striatal location (Redgrave et al., 1999). In the oculomotor loop, for example, the activity in striatum at a particular retinotopic coordinate directly encodes the relative probability of a saccade being made to that coordinate (Hikosaka et al., 2000).

(9) SNr and EP maintain a tonic inhibitory output over their target structures in thalamus and brainstem. Selection of an action is then encoded by suppression of the appropriate neural population (output "channel") in SNr/EP, resulting in selective disinhibition of basal ganglia output targets (Chevalier et al., 1985; Chevalier and Deniau, 1990).

(10) Within the BG, there are many candidate mechanisms with the capacity for selection. At the internucleus level, we included in our model the circuit comprising the output nuclei (represented by the SNr), onto which focused inhibition from D₁-dominant striatal projection neurons is superimposed on diffuse excitation from the STN. This basic scheme was that proposed by Mink and Thach (1993) and Nambu et al. (2002) and developed quantitatively by Gurney et al. (2001a,b). This circuit represents an off-center, on-surround feedforward network (its hypothetical operation is illustrated in Fig. 1). In principle, sufficient focused inhibition from striatum (from the channel with the most salient input) will inhibit the corresponding output channel, resulting in the disinhibition of the targets of that channel. At the same time, diffuse excitation from STN could increase the tonic inhibitory output from the other output channels, enhancing the contrast between selected and nonselected channels. We therefore refer to the circuit formed by the D₁-dominant striatum, STN, and SNr as the "selection pathway."

(11) However, to achieve the correct balance between excitation and inhibition in the selection pathway, excitation from STN needs to be regulated by inhibition from GP. Our simulation and analysis of a population-level BG model showed that the circuit comprising D₂-dominant striatal projection neurons, STN, and GP provided the required amount of inhibition to enable signal selection and switching (Gurney et al., 2001b). We thus termed this circuit the "control pathway." This alternative functional architecture (selection and control pathways) contrasts with the then prevailing qualitative model of basal ganglia anatomy that dissects the system into "direct" and "indirect" pathways (Albin et al., 1989).

(12) In the present study, we are primarily interested in the tonic function of dopamine, associated with pacemaker-like firing in substantia nigra pars compacta (SNc) (Grace, 2002). Thus, we do not include mechanisms for modeling phasic dopamine release by SNc (or the ventral tegmental area) (Schultz, 1998; Dommett et al., 2005), which appears necessary for learning in reinforcement-based contexts. However, because the control of tonic SNc firing by intrinsic basal ganglia connectivity is comparatively poorly understood, we omitted SNc input as a substantive part of the model and incorporate levels of tonic dopamine as a control variable. The effects of dopamine on synaptic transmission are modeled for striatum, STN, and GP. We specifically include effects in STN and GP because of the noted changes in coupling between these nuclei after experimental manipulations that alter tonic dopamine (Magill et al., 2001; Baufreton et al., 2005).

The selection/control BG architecture just described is consistent across the so-called "motor" and "association" macroscopic channels of the BG (Joel and Weiner, 2000) but not, on current evidence, with the "limbic" channel that receives input from prefrontal cortex and outputs from ventral pallidum, without an obvious selection/control pathway distinction. Nonetheless, parts of what may be considered limbic BG, including the nucleus accumbens, have a circuitry remarkably similar to that depicted in Figure 1 (Maurice et al., 1997). Thus, although we specifically restrict our attention here to the selection of motor-based actions, it is conceivable that the same circuitry could cover the spectrum of selection functions, from individual limb movements to cognitive processes (Houk and Wise, 1995; Middleton and Strick, 2000).

Model architecture
Each BG component of the functional anatomy in Figure 1 defines a neuron population in the model (five populations in all). For all simulations reported in this paper, we used N = 3 channels, each with n = 64 neurons per channel, making a total of 192 neurons per population. We denote these populations by the sets _s1 (D₁-dominant striatum), _s2 (D₂-dominant striatum), _T (STN), _G (GP), and _O (SNr). The notation i _O thus means a neuron i in the SNr population, and membership is similarly defined for the other populations.

The channels were assumed to form adjacent neural subpopulations within each of these larger populations. The comparatively small number of neurons in this model was chosen for computational tractability, but this prevented a veridical representation of neuronal population ratios. Therefore, to model the funneling of connections from the relatively large striatal neuron population to its smaller efferent nuclei (Oorschot, 1996), we used large striatal connection weights to emulate the effect of the massive striatal convergence on its target neurons.

Internucleus and intranucleus connectivity was defined stochastically. For most internucleus connections, the potential target neurons were sampled from within the same channel as the afferent neuron. For such connectivity, each neuron in an afferent structure contacted every neuron in the same channel of a target structure with probability ^c, given in Table 1. The exception was STN output, which sampled potential target neurons across all channels in GP and SNr, modeling the relatively diffuse nature of its axon collateralization (Parent and Hazrati, 1995). To ensure comparable mean numbers of efferents across all nuclei, the probability of each STN neuron making contact with one in its target structures was set to ^b = ^c/N.

View this table: [in this window] [in a new window]   Table 1. General parameter values

 
In some experiments, we modeled SNr and GP local axon collaterals (Deniau et al., 1982; Kita and Kitai, 1991). These also sampled potential target neurons across all channels within the same structure, including the parent channel of the neuron in question. This is consistent with the dense local axon collateralization in close proximity to the parent neuron within these structures (Mailly et al., 2003; Sadek et al., 2005): because we modeled channels as adjacent subpopulations, we expected that the local collaterals would extend across all channels. Once again, the probability of making a connection with a target in this diffuse projection scheme was ^b. Each internucleus and intranucleus connection was assigned a transmission delay taken from published data (see Table 3).

Input to the model
Cortical input was emulated by generating spike trains with a given frequency spikes/s. The temporal distribution of cortical spikes has been shown to approximate that derived from a Poisson process in numerous experimental preparations and behavioral contexts (Dayan and Abbot, 2001). For simulations in which normal tonic firing of corticostriatal pyramidal neurons was required, interspike intervals (ISIs) were therefore generated from an exponential distribution with µ = 1/ s, generating a spike train described by a Poisson process (within the temporal quantization determined by the simulation time step t).

For some experiment replications, we required simulations of cortical EEG slow-wave activity. The <1 Hz EEG wave corresponds to alternating periods of silence and an active phase of synchronized tonic firing of cortical pyramidal cells in both anesthetized and sleeping animals (Steriade et al., 1993; Amzica and Steriade, 1995; Steriade et al., 2001). During slow-wave sleep, these cells fire at an overall mean rate of 12 spikes/s (Steriade et al., 2001). This implies that the active phase had a mean rate of at least 24 spikes/s during slow-wave sleep. Given that urethane and ketamine anesthesia result in a greater amplitude <1 Hz cortical EEG oscillation than slow-wave sleep, we expect the corresponding active phase firing rate to be higher (Magill et al., 2001; Steriade et al., 2001). Moreover, action potentials fired during the active phase do not follow a Poisson process (Stern et al., 1997).

To emulate these data, cortical spike trains for all channels were constructed that alternated 0.5 s periods of silence with 0.5 s periods of regularly spaced spikes. For each spike train and for each firing period in that train, a spike frequency f_s was selected from a Gaussian distribution with mean ± SD of 32 ± 6.7 spikes/s. Each spike was then "jittered" in time by 1/f, where f was drawn from a Gaussian distribution with mean ± SD of 0 ± 2.5f_s (with any resulting coincident spikes removed). The resulting cortical input spike trains were not identical and yet non-Poisson (Stern et al., 1997).

Neuron models
We used leaky integrate-and-fire model neurons with additional mechanisms to capture, phenomenologically, key properties of each neuron species. The properties common to all of the neuron models are detailed here: additional components required for each neuron species are described in the Appendix. In general, the neuron membrane potential V^m is given by

where R is the neuron input resistance, I(t) is a sum of various current sources (including synaptic sources), and ^m is the membrane time constant.

A spike is generated when the membrane potential reaches some threshold. To model the repolarization of the neuron, V^m is immediately set to the resting potential V^r after a spike event. With no driving current in Equation 1, V^r = 0; this shift of the resting potential to a convenient value has no functional consequences for the model. We force an absolute refractory period by stopping the solution of Equation 1 for a period ^abs after spike generation.

To model the postsynaptic current (PSC) contribution made to the membrane potential V^m by a single synaptic event, I(t) will contain a contribution I^s(t), defined by a suitably weighted kernel function F(·):

Here, t^f is the time of the presynaptic spike generation, and t^d is the transmission time delay.

The weighting in Equation 2 is conveniently divided into two factors w and Î^s. The latter is the step change in synaptic current caused by a presynaptic spike and therefore could be deduced from the amplitude of the appropriate recorded PSC. However, most synaptic events are recorded as postsynaptic potentials (PSPs) in the whole-cell membrane potential. We therefore fix Î^s so that a single postsynaptic spike event (with w = 1) produces a peak PSP, constrained by data from the literature (see Table 1). The term w is a dimensionless weighting factor that allows us to use each model synapse to represent multiple physical synapses.

The kernel function F(·) is a step followed by an exponential decay (Gerstner, 1999), which is a simple fit to recorded PSCs:

where H(x) is the Heaviside step function [H(x) = 0 for x < 0 and H(x) = 1 otherwise], and ^s is the synaptic time constant (a measure of the characteristic duration of the PSC).

We model the three functionally predominant types of synaptic currents in basal ganglia: fast EPSCs and IPSCs attributable to AMPA and GABA_A receptors and slow EPSCs attributable to NMDA receptors. [For simplicity, we consider only direct synaptic actions of glutamate and GABA and omit GABA_B receptors in the present model: postsynaptic activation of these receptors induces a intracellular signaling cascade rather than a synaptic PSC (Bormann, 1988).] These give rise to synaptic currents I^ampa(t), I^gabaa(t), and I^nmda(t), respectively, with associated peak PSPs of V^ampa, V^gabaa, and V^nmda. For each of the synaptic types modeled, we have separate synaptic time constants (^ampa, ^gabaa, and ^nmda). (This is a simplistic model of the NMDA receptor-evoked current, omitting both its short rise time, as is often done, and its voltage-dependent threshold for activation.) The total synaptic current is then found by summing over all afferents and all postsynaptic events up to the current time t. For a generic synaptic type, with time constant ^s, this is given by

where is the set of afferents to the target neuron, _j is the set of firing times t^f for afferent neuron j, and t is the transmission delay between j and the target neuron.

Synaptic efficacy is modulated by dopamine within the striatum, STN, and GP (for details of individual models, see the Appendix). We parameterize tonic dopamine concentration at D₁ and D₂ receptors by _D1 and _D2, respectively (with 0 _D1, _D2 1). By separately parameterizing contributions to D₁ and D₂ receptors, we were able to mimic the effects of receptor-specific agonists and antagonists if required. Under normal conditions, _D1 = _D2.

An important determinant of synaptic efficacy is the relative location of the synapses on the postsynaptic membrane. Inhibitory synapses on proximal dendrites or on the soma are postulated to provide a strong shunting or "vetoing" inhibition effect that gates input from more distal dendritic sources (Blomfield, 1974; Shepherd and Brayton, 1987; Segev, 1995; Ulrich, 2003). That is, the effect of activating postsynaptic receptors at proximal or somatic synapses is more divisive than subtractive, thereby attenuating distally originating PSPs. The pattern of synapse location within the basal ganglia suggests that this mechanism may be important within STN, GP, and SNr.

To model the shunting effect of proximal and somatic inhibition, we adopt a phenomenological description, or "pseudocompartmental" scheme, similar to that used in our previous modeling study on oscillations in GP and STN (Humphries and Gurney, 2001). When shunting inhibition is present, our model contains three pseudocompartments, corresponding to distal and proximal dendrites and the soma. Synaptic input from each GABAergic afferent j is randomly assigned to one of these compartments with probabilities P(s), P(p), and P(d) (for somatic, proximal, and distal compartments, respectively). These probabilities are based on estimates from the literature for the relative proportions of GABAergic contacts within each zone of the target neuron species. The result is that, for any neuron, each of its GABAergic afferents are assigned to either a somatic , a proximal , or a distal compartment set.

Excitatory inputs within basal ganglia are supplied by STN neurons whose axon terminals mainly contact the distal dendrites of their target neurons in SNr and GP (Bevan et al., 1994; Shink et al., 1996). Excitatory inputs from cortex to STN and striatal projection neurons also mainly contact their distal dendrites (Chang et al., 1983; Bevan et al., 1995; Gerfen and Wilson, 1996). Thus, all excitatory afferents are effectively assigned to for simplicity.

The distal synaptic current I_i^D going into the proximal compartment on neuron i is computed additively, as follows:

where each of the terms on the right side is given by an expression of the form of Equation 4 and where notation of j _i limits the scope of afferent set in Equation 4 to just those afferent GABAergic synapses that fall in the distal compartment.

One potential problem with a current source-based scheme and additive inhibition is that massive inhibition can hyperpolarize the membrane to biologically unobtainable values. This will also prevent subsequent excitation from working realistically. We therefore limit V^m from below by a threshold V^lim, which mimics the effect of the reversal potential E_Cl for chloride ions on GABA-elicited PSPs. E_Cl is typically –70 mV, with measured resting potentials of –50 mV. Therefore, under the shift of resting potential by 50 mV to give V^r = 0, V^lim was set to –20 mV.

We define the effects of shunting inhibition on the distal input in two variables h and h, for the contributions of proximal and somatic inhibition to the ith neuron, respectively. These take values in the interval [0,1], where a value of 0 indicates a complete veto of all input from the compartments more distant from the soma (a full definition is given in the Appendix). Shunting inhibition simultaneously forces the membrane potential toward the reversal potential associated with GABA_A receptors, E_Cl (which, under the shift of resting potential V^r to 0, is equivalent to V^lim in the model). To model this, we introduce a constant current I, which drives the membrane potential to V^lim if fully activated and which is gated by a variable Q, which is dependent on h and h. Again, Q takes a value in the interval [0,1], where a value of 1 indicates that I is fully activated (see the Appendix).

Every neuron species has some form of constant input current, I^spon, which models the effect that intrinsic membrane currents have on spontaneous activity. Furthermore, all neuron models contain a noise term, which subsumes noise produced by several sources, including graded vesicle release and failure (Manwani and Koch, 1999). Noise is added at every time step as a voltage deflection of V^m sampled from a Gaussian distribution with mean ± SD of 0 ± 0.3 mV. Finally, it is occasionally useful to incorporate currents I*(t) that model specific time-varying phenomena. Specifically, we incorporate a piecewise approximation [labeled I*(t) = I^Ca] to the current activation sequence that underlies rebound burst firing in STN neurons (Beurrier et al., 1999) (we define this in the Appendix). Incorporating all contributions defined above, the membrane equation becomes (dropping index i for clarity)

In the special case in which the neuron class does not have proximal or somatic inhibition, we have no need for the pseudocompartmental scheme and write

where I^s is the synaptic input. The basic membrane Equation 6 is adapted for each neuron species according to specific models of dopaminergic modulation and additional properties of those species. The adaptations are described in detail in the Appendix.

Modeling strategy
Models-as-animals. We briefly address here the design of the numerical simulations reported in Results. The statistical properties of experimental data are dependent on both the number N_D of data points collected and the manner in which the data points are collected. Experimental studies routinely pool cell recordings taken from different animals, and, if we are to compare the model and the data in a principled way, we need to replicate this data-gathering process. Our models are structurally stochastic, and so we advocate a "models-as-animals" approach in which each instantiation of the model is considered to represent a single animal. Then, for each experimental condition, we construct the same number of models N_A as there were animals and collect n_D observations (i.e., cells) from each neural structure in that condition, so that N_An_D is as close to N_D as possible (this will not necessarily be identical because automated postprocessing relies on uniform sampling). This constitutes a single "virtual experiment." We then repeat this many times to determine the spread of data that the model predicts for each measured experimental variable.

In this way, we aim to generate results whose statistical properties are comparable with their experimental counterparts, thereby enabling statistically valid claims about the "fit" between the model and the data. This approach should be contrasted to that more commonly used method in which a single instantiation of a model is simulated, the appropriate model parameters altered to mimic the experimental manipulations, and all the outputs of all the cells in the model analyzed (additional discussion of these issues has been given by Humphries and Gurney, 2006).

Determination of parameter values. Our BG model was explicitly intended as an in vivo network. Every parameter specific to a neuron species (e.g., R, ^m) had an experimentally determined value, taken from in vivo studies in which those were available. The values and their literature sources are given in Table 2. General neuron parameters took biologically plausible values based on those presented previously (Dayan and Abbot, 2001) and are given in Table 1.

View this table: [in this window] [in a new window]   Table 2. Population-specific parameters

 
Moreover, we included potentially crucial biological detail often omitted by modelers. Any model that intends to study oscillatory or correlation phenomena must incorporate interneuronal transmission delays, because these critically determine the timing of spikes. In our model, every connection had a transmission delay with experimentally determined values. The distribution of synapses from each connection to the somatic, proximal, and distal pseudocompartments of a neuron species was also taken from the experimental literature. The delay values and synapse distributions, along with their literature sources, are given in Table 3.

View this table: [in this window] [in a new window]   Table 3. Connection parameters

 
Furthermore, to avoid spurious correlations and oscillations between neurons, critical parameter values within a neuron population were assigned stochastically. The resistance R and membrane time constant values ^m for each member of a neuron population were sampled from a Gaussian distribution, with a mean given by the appropriate values in Table 2 and SD that was 10% of that mean value (giving a conservative, biologically plausible, coefficient of variation of 0.1). This ensured a sufficient range of all these parameters values so that neurons within the same structure did not fire identically for a given input (the burst-firing current parameters for STN neurons were similarly distributed).

The only model parameters that could not be assigned directly from experimental data (relative weighting of shunting inhibition and overall network connection density) were set at the outset and then not changed for any of the simulations reported here. Tonic input currents for STN, GP, and SNr neurons were set to give the approximate tonic firing rates for those nuclei that matched experimentally measured rates. This provided a baseline version of the model (see Results). The only parameters changed thereafter were those required to emulate the conditions of the experimental study we were simulating and are fully described in Results.

Data analysis
The firing rates of individual neuron outputs were computed by a moving-window average, using an -function window of 10 ms width computed every 1 ms (Nawrot et al., 1999; Dayan and Abbot, 2001). Autocorrelograms were computed with 10 ms bins, following standard methods (Abeles, 1982; Dayan and Abbot, 2001). Power spectra were computed directly from the spike trains using the multi-taper spectrum method, which is particularly suitable for analyzing discrete point time series such as spike trains (Jarvis and Mitra, 2001). The Matlab (MathWorks, Natick, MA) toolbox Chronux, which implements multi-taper methods, is available from http://chronux.org/. Cortical pseudo-EEG oscillations were analyzed using the Lomb–Scargle periodogram (Scargle, 1982). All analysis functions, other than the Chronux toolbox, were custom written in Matlab.

Implementation
All simulations were implemented using Matlab and custom MEX-file components in C for efficiency. They were run on a Compusys cluster of dual Opteron servers. The code is available on request or from our website www.abrg.group.shef.ac.uk. All V^m equations were solved exactly (with the weak assumption that I^Ca is constant over a small time interval) and simulated in discrete time to catch all afferent spike arrival and neuron firing events. The simulation time step was t = 0.1 ms throughout.

	Results

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Tonic firing can be fitted from in vivo recordings of alert rats
The first step in using the model was to "calibrate" it so that its quiescent state gave realistic tonic firing rates. However, given the complexity of the model, it is not clear a priori if it is possible to fit the tonic rates to target data. The ability to do so therefore constitutes a result from the model.

The hypothesis concerning the role of the BG in action selection in the behaving animal is central to our work. We therefore chose to adopt a resting state consistent with data from non-anesthetized animals. Thus, target firing rates were taken from recordings in awake, resting rats of cells in STN (Kreiss et al., 1997), GP (Urbain et al., 2000), and SNr (Windels and Kiyatkin, 2004).

The calibration process consisted of varying the I^spon currents in each of STN, GP, and SNr while applying background cortical input consisting of low-frequency, unpatterned (Poisson series) spikes with a mean firing rate of = 3 spikes/s (Bauswein et al., 1989). Such a low rate is insufficient to cause reliable firing in striatum, consistent with numerous striatal recordings (Mahon et al., 2001, 2003). During this process, each simulation was run for 10 s. Mean firing rates were computed over the last 9 s to avoid the confounding effects of initial transients in model neuron output.

To compare the simulation output with experimentally determined quantities, we "recorded" five model neurons in each of STN, GP, and SNr, from each of 11 instantiations of the model that represented 11 "virtual animals" (see Materials and Methods). [The original experimental papers reporting the mean firing rates (Kreiss et al., 1997; Urbain et al., 2000; Windels and Kiyatkin, 2004) do not state the number of animals used, and so we assumed five cells from each animal, based on typical experimental protocols.] Each batch of 11 model simulations constituted a "virtual experiment," which was repeated 50 times to get a measure of the spread of results we could expect using the given (real and simulated) protocol. Dopamine levels were set to _D1 = _D2 = 0.3 to simulate normal tonic dopamine concentration. These values were used throughout, except in simulations of experimental conditions with altered dopamine (see the following sections).

Figure 2A shows the mean firing rates for the model that included collaterals within SNr and GP, together with the mean rates derived from in vivo recordings. To gain a quantitative measure of how well the model fits the data, we ran a two-tailed independent-groups t test between the experimental data and each virtual experiment (50 t tests in all). For each of STN, GP, and SNr, the percentage of virtual experiments yielding mean tonic firing rates that were statistically indistinguishable (at p = 0.05) from the experimental data were 84, 92, and 40%, respectively.

View larger version (7K): [in this window] [in a new window]   Figure 2. Mean tonic firing rates from STN, GP, and SNr neurons in the awake resting state of the simulated and real basal ganglia. The box-and-whisker plots summarize the distribution of mean firing rates across the 50 virtual experiments. The mean firing rates from the biological experimental data (see Results) are shown by the open circles, with error bars (thick lines) showing the 95% confidence interval for each mean value. A, Model with collaterals in SNr and GP. B, Model with no collaterals in SNr and GP. Both forms of the model were able to reliably capture the respective mean firing rate properties of these nuclei in the awake rat.

 
Figure 2B shows the corresponding results for a model with no local axon collaterals within SNr and GP (part of the program of work detailed in a later section). The percentages of nonsignificant t tests in this case was 60, 100, and 86%. With or without collaterals, the mean tonic firing rates of the model nuclei are mostly indistinguishable from the experimental data.

For the full model, with GP and SNr collaterals, the resulting values for the I^spon currents were as follows: for STN, I = 11 x 10^–10 A; for GP, I = 3.8 x 10^–10 A; and for SNr, I = 3.9 x 10^–10 A. These values were used for all of the following simulations, except when noted below. Without collaterals in the GP and SNr, the resulting values for the I^spon currents were as follows: for STN, I = 9 x 10^–10 A; for GP, I = 3 x 10^–10 A; and for SNr, I = 3.4 x 10^–10 A. These values were used for one batch of simulations detailed in the section below on slow-wave oscillations.

To minimize initial transients for future simulations, a complete snapshot of the state of the model, including membrane potential values and synaptic current totals, was recorded at the end of one of the simulations described here and was used as the initial conditions for all additional simulations.

Input signals are appropriately selected and switched between
We first established whether the model retained the capacity to support the functional hypothesis of action selection. To assess the selection and switching performance quantitatively, we followed the procedures used in our previous work (Gurney et al., 2001b), so that a direct comparison could be made with the population-level model described there. Thus, all simulations in this section used the following stimulus protocol. An input to channel 1 with mean firing rate (1) was initiated at t = 1 s, followed at t = 2.5 s by an input to channel 2 with rate (2); the simulation terminated at t = 5 s. This sequence gives three time intervals: I₁ = [0, 1), I₂ = [1, 2.5), and I₃ = [2.5, 5]. The magnitude of each input represents the salience of the requested motor program: hence, the sequence used here tests the ability of the BG to select an action (at t = 1) and then to switch to a more salient action (at t = 2.5) but only if it is appropriate to do so.

The model was run with combinations of (1), (2), drawn from the range [4, 40] spikes/s with step-size 4 spikes/s, giving 10 values of each firing rate and 100 experiments in all. The set of 100 experiments was repeated with different levels of simulated tonic dopamine: _D1 = _D2 = 0, _D1 = _D2 = 0.3, and _D1 = _D2 = 1, corresponding to "low," "normal," and "high" levels of dopamine, respectively. However, because there was no specific experimental data with which to compare the model, only a single model was run under each condition. Selection was deemed to have taken place if the mean firing rate of an SNr channel was less than some threshold ^s = 5 spikes/s, which signaled sufficient disinhibition (Chevalier and Deniau, 1990).

Figure 3 shows SNr output from three example simulations, one at each dopamine level, from a model that received inputs (1) = 20 spikes/s and (2) = 40 spikes/s. The output shown is the mean instantaneous firing rate across all neurons in a particular SNr channel. The first example is from the normal dopamine condition (Fig. 3A). During I₁, all three channels show tonic firing in the quiescent state. After the onset of a salient input to channel 1 of the BG during I₂, the mean output of SNr channel 1 falls below ^s, and the mean output of the other channels increases above the original tonic level. Subsequently, after the onset of a more salient input to channel 2 of the BG during I₃, the mean output of SNr channel 2 falls below ^s, whereas, at the same time, the mean output of SNr channel 1 returns to a level above ^s. Thus, the mean SNr outputs correctly signaled the selection of the most salient action and then correctly switched to the selection of a more salient action once it had become available. The histograms on the right of Figure 3 show the distributions of individual SNr neuron mean firing rates within channels 1 and 2. The changes in the distribution shape reflect the selection and switching process, yet they also show that not all neurons in a channel responded in the same manner, and a wide range of outputs were maintained although the mean output of the channel was signaling the correct response. This result highlights the difficulty of pinpointing the channel (loop) membership of an individual neuron from its firing properties alone.

View larger version (25K): [in this window] [in a new window]   Figure 3. Examples of selection and switching under different levels of simulated dopamine. In each case, the large panel on the left shows the instantaneous mean firing rate for each SNr channel, obtained by averaging across every SNr neuron in the channel. Below this, the rectangular waveforms indicate mean firing rates of cortical inputs driving the corresponding channels in striatum and STN. The small panels on the right show histograms of individual mean firing rates of neurons in channel 1 (Ch1) and channel 2 (Ch2) during each time interval. These histograms thus show the affect on the whole channel population of each change in input. A, Normal levels of simulated tonic dopamine (D1 = D2 = 0.3). B, Low levels of simulated tonic dopamine (D1 = D2 = 0). C, High levels of simulated tonic dopamine (D1 = D2 = 0.8).

 
The second example, in Figure 3B, shows how low simulated dopamine typically results in no selection. Neither of the salient inputs to the BG is sufficient to cause SNr channel output to fall below ^s. The distributions of the individual neuron firing rates change little in each time interval compared with the normal dopamine condition (in Fig. 3A), showing that the whole channel population did not have a clear response to the input. The third example, in Figure 3C, shows how high simulated dopamine can result in the simultaneous (here dual) selection of salient inputs. Both SNr channels 1 and 2 are selected in interval I₃, because the output of SNr channel 1 did not rise above ^s after the onset of salient input to channel 2 of the BG. This is reflected in the highly skewed distribution of individual mean firing rates for neurons in SNr channel 1 during I₃. The outcome of all three example simulations closely matches the corresponding behavior of our previous population-level BG models (Gurney et al., 2001b) and corroborates the hypothetical selection mechanism discussed in Materials and Methods.

The outcome of each of the 100 experiments for each model was classified according to the mean SNr firing rate for channels 1 and 2 in intervals I₂ and I₃. Bearing in mind that selection of a channel is deemed to have occurred if the mean SNr firing rate of that channel is below the selection threshold ^s, the following possibilities were recognized: (1) no selection: neither channel is ever selected in any interval (as in the example of Fig. 3B); (2) selection: a single channel is selected during the competition; that is, either channel 1 is selected in both intervals I₂ and I₃ and channel 2 is never selected, or channel 1 is never selected and channel 2 is selected in I₃; (3) switching: channel 1 is selected in I₂ and deselected in I₃ with subsequent selection of channel 2 in I₃ (as in the example of Fig. 3A); (4) dual selection: both channels are selected in I₃ (as in the example of Fig. 3C); and (5) interference: all other cases.

A summary of the SNr output categorization across all tested input pairs for a single model under normal dopamine is shown in Figure 4A. The majority of input pairs resulted in successful resolution of the competition between the two channels, that is, single selection or switching. With smaller magnitude inputs, there was no selection (bottom left of grid). This behavior, corresponding to a filtering of low salience inputs, is primarily attributable to the existence of a down state of the striatal neurons, as modeled by the I^down current (see the Appendix). The cutoff of the selection filter is 16 Hz. These results closely match those of our previous population-level models (Gurney et al., 2001b, 2004b).

View larger version (40K): [in this window] [in a new window]   Figure 4. Selection and switching capabilities of the basal ganglia model. Each simulation was assigned to a category based on its combination of mean SNr outputs on channel 1 and channel 2 in response to the inputs (see Results). A, Normal levels of tonic dopamine (DA) (D1 = D2 = 0.3). B, Low levels of simulated dopamine (D1 = D2 = 0.0). C, High levels of simulated dopamine (D1 = D2 = 0.8). D–F, Idealized templates for selection and switching under normal, low, and high dopamine levels, respectively (see Results). G, Box-and-whisker plots of the matches between a batch of 50 models and corresponding idealized selection templates, for each level of simulated dopamine.

 
We also tested the effects of both dopamine depletion and excess dopamine on the selection and switching of input signals. To simulate dopamine depletion, we set _D1 = _D2 = 0 and ran the same model used with normal dopamine levels. With low dopamine levels, no signal selection could occur, and no switching was observed (Fig. 4B). Again, these results replicate those of our population-level model (Gurney et al., 2001b). The implication of these results is that, after dopamine depletion, actions of low urgency are not selected by the basal ganglia.

To simulate excess dopamine, we set _D1 = _D2 = 0.8 and ran the same model again. Switching now rarely occurs: competition between salient inputs now most often results in dual-channel selection (Fig. 4C). This also replicates results from the population-level model (Gurney et al., 2001b).

To quantify how well the model performs in these experiments, we defined templates of the idealized outcomes for the selection competition under different dopaminergic conditions (Fig. 4D–F). These templates define different selection/switching regimens for a single-action selection mechanism that are under the control of dopamine. With normal tonic dopamine, the outcome for such an idealized selection mechanism is as follows: no selection if both inputs are less than the selection filter cutoff; selection of a single channel if its input is above the filter cutoff and the other channel has a filtered input; or, if both channel inputs are above the filter cutoff, switching to channel 2 after the onset of its input only if it has a higher salience than the input to channel 1. The result is the idealized template in Figure 4D. With all tonic dopamine depleted, the idealized outcome is no response to any salient input in the tested range (Fig. 4E). With high levels of tonic dopamine, the idealized outcome is either no selection if both inputs are less than the selection filter cutoff or dual selection otherwise (Fig. 4F): if both actions are salient, then no clear decision is made between them. [The rationale for these regimens and their construction has been discussed previously (Gurney et al., 2004b).]

The performance of a model over the set of 100 experiments with a particular level of dopamine may now be quantified by calculating its percentage match with the template (via the fraction of matching outcomes). This was done for 50 instantiated models (i.e., 50 sets of 100 experiments for each model), for each level of dopamine (low, normal, and high). There is a good match with each of the idealized templates for every tested model (Fig. 4G), showing that the model consistently performs appropriate selection and switching for the simulated level of tonic dopamine.

Cortical slow-wave activity and the STN–GP loop
We now validate the model by replicating detailed experimental data and derive novel predictions. A striking result from recent studies of the STN–GP loop is the demonstration of the role of dopamine in controlling the relative independence of the firing patterns of these nuclei (Magill et al., 2001) (see Fig. 7C,D). Under urethane anesthesia, most STN neurons show burst firing with an interburst frequency of 1 Hz. This low-frequency oscillation (LFO) is highly correlated with the cortical EEG slow wave, which also has a frequency of 1 Hz. GP neurons are not correlated with the cortical EEG and display no LFOs. Cortical ablation stops the LFOs in the STN but does not alter GP behavior. Thus, STN activity appears to be driven by cortex, but, contrary to intuition, the strong oscillation and changes in firing pattern are not transmitted to the GP, despite the direct excitatory pathway from STN to GP.

After dopamine depletion (by 6-OHDA lesion of the SNc in otherwise intact animals), most GP neurons do show LFOs that are correlated with the cortical EEG slow wave. This suggests that some action of dopamine prevents cortical EEG slow-wave activity from being transmitted to the GP via the STN in the intact animal. (We particularly sought to investigate this phenomenon and, therefore, included in the model dopaminergic modulation of inputs to the STN and GP, as described in the Appendix.)

In addition, Magill et al. (2001) found that combining the ablation of cortex and the lesion of SNc left residual LFO activity in both STN and GP despite the absence of cortical input. This suggests the possible presence of a natural pacemaker circuit formed by STN and GP (Plenz and Kitai, 1999). The underlying mechanism for this pacemaker could take the form of that proposed by Humphries and Gurney (2001), which relies on GP-induced rebound bursting in STN; this proposal was therefore also investigated with the model.

Replication of these results also requires consideration of the role of anesthesia. The "baseline" state of the model described in the previous sections was constrained by data from awake, unrestrained rats. If the model is to explain the LFO phenomena, it must also match firing rates under urethane anesthesia. Recent work has shown that urethane acts to increase the efficacy of GABA synapses but decrease the efficacy of glutamate synapses (Hara and Harris, 2002). This was modeled by multiplying all inhibitory connection weights by 1.5 and all excitatory weights by 0.65, matching the effect magnitude reported previously (Hara and Harris, 2002). We also altered the spontaneous current in STN, I, to 5 x 10^–10 A. These parameter values were determined by matching the mean firing rates of STN and GP cells in the control condition only (see below) to those reported previously (Magill et al., 2001). Alteration to I was necessary to achieve this match, because it was clear that urethane anesthesia alters the spontaneous activity of STN cells compared with awake animals [compare the experimental data from the awake rats (Fig. 2) with data from the anesthetized rats (see Fig. 7A, control condition)].

Four sets of simulations were run, emulating the design and manipulations of the four experimental conditions (Magill et al., 2001):

(1) The control condition included an intact cortex and basal ganglia: slow-wave input (see Materials and Methods) and _D1= _D2 = 0.3, representing normal tonic levels of dopamine. This condition had N_A = 7 animals, sampling n_D = 4 STN cells from each animal (giving a total of 28 STN cells), and n_D = 3 GP cells from each animal (giving a total of 21 GP cells).

(2) The ablation condition included an ablated cortex and intact basal ganglia: corticostriatal weights set to 0 and dopamine parameters _D1= _D2 = 0.3. This condition had N_A = 7 animals, sampling n_D = 2 STN cells (giving a total of 14 STN cells), and n_D = 4 GP cells (giving a total of 28 GP cells).

(3) The lesion condition included an intact cortex and SNc lesion: slow-wave input and with loss of tonic dopamine _D1= _D2 = 0. This condition had N_A = 6 animals, sampling n_D = 6 STN cells (giving a total of 36 STN cells), and n_D = 5 GP cells (giving a total of 30 GP cells).

(4) The lesion-and-ablation condition included an ablated cortex and SNc lesion: corticostriatal weights set to 0 and _D1= _D2 = 0. This condition had N_A = 6 animals, sampling n_D = 3 STN cells (giving a total of 18 STN cells), and n_D = 3 GP cells (giving a total of 18 GP cells).

All simulations were run for t = 10 s. Our analyses followed those performed by Magill et al. (2001) on their data. However, because we had no separate EEG per se, we created a pseudo-EEG for the cortical activity to enable computation of the spike-triggered averages. Because the cortical EEG corresponds to the synchronous activity of a large population of cortical neurons (Stern et al., 1997; Steriade et al., 2001), the pseudo-EEG for each simulation was constructed by finding the firing rate of each cortical spike train using a moving-window average (see Materials and Methods) and then averaging these firing rate reconstructions over all cortical inputs. Spike-triggered averages of the pseudo-EEG signal were then computed using a window of ±1 s from each spike.

The model replicates individual spike train properties from the experimental data
The two most striking, non-intuitive results of Magill et al. (2001) were the complete decoupling of STN and GP activity in the control condition and the presence of LFO activity in both nuclei in the combined lesion-and-ablation condition. In Figures 5 and 6, we show that the model can replicate the detailed firing properties of individual STN and GP neurons reported by Magill et al. (2001)