Abstract
The basal ganglia (BG) play an important role in motor control, reinforcement learning, and perceptual decision making. Modeling and experimental evidence suggest that, in a speed–accuracy tradeoff, the corticostriatal pathway can adaptively adjust a decision threshold (the amount of information needed to make a choice). In this study, we go beyond the focus of previous works on the direct and hyperdirect pathways to examine the contribution of the indirect pathway of the BG system to decision making in a biophysically based spiking network model. We find that the mechanism of adjusting the decision threshold by plasticity of the corticostriatal connections is effective, provided that the indirect pathway counterbalances the direct pathway in their projections to the output nucleus. Furthermore, in our model, changes within basal ganglia connections similar to those that arise in parkinsonism give rise to strong beta oscillations. Specifically, beta oscillations are produced by an abnormal enhancement of the interactions between the subthalamic nucleus (STN) and the external segment of globus pallidus (GPe) in the indirect pathway, with an oscillation frequency that depends on the excitatory cortical input to the STN and the inhibitory input to the GPe from the striatum. In a parkinsonian state characterized by pronounced beta oscillations, the mean reaction time and range of threshold variation (a measure of behavioral flexibility) are significantly reduced compared with the normal state. Our work thus reveals a specific circuit mechanism for impairments of perceptual decision making associated with Parkinson's disease.
Introduction
The basal ganglia (BG), working together with the cortex, are crucial for action selection and decision making (Graybiel, 1995; Mink, 1996; Ding and Gold, 2013). In a decision process, a motor response could be triggered when accumulated information about alternative choices reaches a threshold level. A computational study found that the corticostriatal connection can adjust the threshold over a wide range (Lo and Wang, 2006). The proposal, which has gained empirical support from human imaging experiments (Forstmann et al., 2008), demonstrates that speed–accuracy tradeoff in decision making can be implemented by flexibly adjusting the threshold, and the corticostriatal pathway is modifiable through dopamine-dependent synaptic plasticity (Surmeier et al., 2010). In the proposed mechanism, neural activity in the BG output nucleus, the substantia nigra pars reticulata (SNr), remains unchanged until an abrupt decrease takes place, leading to disinhibition of the superior colliculus (SC) that combines with cortical excitation to trigger a saccadic response. This nonlinearity is important for the corticostriatal pathway to produce a wide range of decision threshold levels.
Recent experiments, however, found gradual ramping of neuronal activity in the caudate (CD) nucleus before a decision threshold is reached in monkeys performing a perceptual decision task (Ding and Gold, 2010). Besides revealing an active involvement of the CD in the perceptual decision process, the observation of ramping activity also calls for generalization of the threshold modulation mechanism suggested previously (Lo and Wang, 2006). With only the direct pathway incorporated in the model, ramping activity at the CD would lead to increasing inhibition of SNr neurons, which would gradually ramp down rather than displaying a sudden drop of activity, raising the question of whether the main conclusions of Lo and Wang (2006) would still hold in that scenario.
The classical model for the BG architecture suggests that normal function requires a balance of the direct pathway, from the striatum to BG output nuclei, with the indirect pathway, which includes a reciprocal loop composed of the subthalamic nucleus (STN) and the external segment of globus pallidus (GPe; Albin et al., 1989; Hikosaka et al., 2000; Calabresi et al., 2014). The main purpose of this work is to investigate the role of the indirect pathway in perceptual decision making. In particular, we generalize the previous threshold detection and modulation mechanism (Lo and Wang, 2006) by developing a decision-making circuit model incorporating the indirect pathway of the BG and the experimentally observed ramping activity at the CD.
The BG are centrally implicated in parkinsonism, one important feature of which is enhanced oscillations in the population neuronal activity in a broadly defined beta band (8–30 Hz; Boraud et al., 2005; Gatev et al., 2006; Hammond et al., 2007; Obeso et al., 2008). The strong interconnections between the STN and GPe highlight the role of this subcircuit as a natural source of rhythmic BG activity (Plenz and Kital, 1999; Gillies et al., 2002; Terman et al., 2002; Holgado et al., 2010; Kumar et al., 2011). Our model suggests a circuit mechanism for the generation of beta oscillations in the STN–GPe subnetwork under experimentally identified parameter modulations, which offers novel support for the role of NMDA in this activity pattern. We investigate the influence of the resulting beta oscillations on perceptual decision-making behaviors and show that the mean reaction time and range of threshold values are significantly reduced.
Materials and Methods
Behavioral task simulation.
We constructed and tested our model focusing on a widely used perceptual decision task, the random-dot visual direction discrimination task (Newsome et al., 1989; Roitman and Shadlen, 2002; Ding and Gold, 2010, 2013). In that experiment, partially coherent moving dots are presented and the subject is required to make a two-alternative choice about the coherent direction (the reaction-time version of the task). In our model simulation, the direction of coherent motion is assumed to be leftward or rightward. Two signals, representing sensory input induced by the leftward and rightward moving dots, are fed into left-preferring and right-preferring cortical modules CxL and CxR, respectively, in a model cortical network (Fig. 1). The firing times of each input neuron are sampled from a Poisson process with a mean rate μ. The mean μ depends linearly on the dot coherence level and is given by the equations μ0 + μA × c for the preferred direction and μ0 − μB × c for the nonpreferred direction, where μ0 (20 Hz) is the baseline input, c (between 0 and 100%) is the coherence level, and μA (60 Hz) and μB (20 Hz) are proportionality factors (Wang, 2002). Cortical outputs project to a model SC with two excitatory projection units, SCL and SCR. The decision time is defined as the time interval between the start of the sensory input to the cortex and the onset of activity in an SC unit, SCe, representing a saccadic burst, where the subscript “e” is either L or R, corresponding to the population that becomes active. The decision threshold is defined as the firing rate of the cortical population Cxe that projects to the active SCe unit at the onset time of SCe activity. A correct trial is defined as a trial in which the model generates a saccadic burst in the same direction as the coherent motion. Model performance is defined as the fraction of correct trials.
Schematic circuit model with three interconnected structures: the cortex, BG, and SC. The two selective excitatory populations in the cortical circuit exhibit ramping activity reflecting integration of sensory input and compete with each other through shared inhibition. Sensory evidence accumulated in the cortical network is sent to the SC (a command center for eye movement) and through the BG. The SC sends feedback signals to the cortex. Inset, The structure of the BG used in our model. Signals to the CD are transmitted to the output nucleus (SNr) through two pathways: the direct pathway (Direct) from the CD directly to the SNr and the indirect pathway (Indirect) from the CD through the GPe, STN to the SNr. Exc, Excitatory; Inh, inhibitory; subscript L, left; subscript R, right. Cxi and SCi are inhibitory populations in the Cx and SC, respectively.
Neuron model.
Neurons in the STN and GPe display postinhibitory rebound (PIR) bursts, i.e., a brief burst of spikes is produced at the end of a relatively long (∼100 ms) period of hyperpolarization (Bevan et al., 2002a; Kita, 2007), which is strongly enhanced in the parkinsonian state (Filion and Tremblay, 1991; Bergman et al., 1994; Boraud et al., 1998; Wichmann and Soares, 2006). We use the integrate-and-fire-or-burst model (Smith et al., 2000) to capture this firing characteristic of STN and GPe neurons. The membrane dynamics is given by the following equation (Eq. 1):
where H(V) is the Heaviside step function and where C = 0.5 nF, gL = 25 nS, VL = −70 mV, Vh = −60 mV, and VT = 120 mV. When the membrane potential reaches a boundary Vb, it is reset to Vr, where Vb = −50 mV and Vr = −55 mV. In this model, rebound burst firing is induced by a low-threshold T-type calcium current. The inactivation variable h of the current satisfies the following equations (Eqs. 2 and 3):
and
where gT = 60 nS, τh− = 20 ms, and τh+ = 100 ms. For neurons in the cortex, SC, and other nuclei in the BG, we use the leaky integrate-and-fire model, which is described by Equation 1 with gT = 0.
The synaptic current Isyn is given by the following equation (Eq. 4):
where the reversal potentials VE = 0 mV and VI = −70 mV. gk is the synaptic efficacy, where the indices k = 1, 2, 3 indicate AMPA, NMDA, and GABAA synapses, respectively. The gating variable s satisfies the following equations (Eqs. 5 and 6):
for AMPA and GABAA receptor-mediated currents and:
for NMDA receptor-mediated current, where tj is the time for the jth spike and α = 0.63. For the decay time constant τ, we use 2 ms for AMPA, 5 ms for GABAA, and 100 ms for NMDA-mediated currents. Isyn describes both the synaptic inputs from other neurons within the circuit and background inputs modeled as Poissonian spike trains representing sources beyond the circuit. We also use a short time delay, 0.2 ms, for synaptic transmission in both normal and parkinsonian states.
To check the possible impact of saturation of CD activity during evidence accumulation (Ding and Gold, 2010), we introduce short-term depression (STD) at corticostriatal synapses (Lovinger et al., 1993; Ding et al., 2008). The STD is implemented by including a factor D that multiplies all terms on the right side of Equation 4 (Hempel et al., 2000), which satisfies the following equation (Eq. 7):
where p = 0.45 and τD = 600 ms.
Network structure.
The full circuit includes three brain areas (Fig. 1): the cortical network (Cx), the BG, and the SC. The Cx performs evidence accumulation, as observed, for example, in the lateral intraparietal area and frontal eye field, and is described in previous work (Wang, 2002). The present model extends an earlier one (Lo and Wang, 2006) by including the indirect pathway of the BG, encompassing a projection from the CD (the eye movement part of the striatum) to the STN–GPe loop as well as projections from both the STN and GPe to the SNr. The SNr is the output nucleus of the BG involved in saccadic eye movements (Basso and Sommer, 2011). The projection directly from the cortex to the STN, the so-called hyperdirect pathway, is simplified to a background input to the STN.
Each BG nucleus includes two populations of neurons selective to left and right moving dots, respectively. Each population in the CD and SNr includes 250 neurons with all-to-all connections. Each population in the STN and GPe includes 2500 neurons, since we will use sparse connectivity within the STN–GPe subcircuit to prevent the population spike pattern from being oversynchronized in the parkinsonian state. Following another model (Kumar et al., 2011), the connection probabilities used are 0.05 for STN→GPe, 0.05 for GPe→GPe, and 0.02 for GPe→STN, reflecting the sparse connections between the STN and GPe and also within the GPe (Kita and Kitai, 1994; Bevan et al., 1997; Sadek et al., 2007; Baufreton et al., 2009). There are no known recurrent connections within the STN (Sato et al., 2000). Hence, such connections are not included. All other projections are all-to-all for the left and right selective populations, respectively. The SC network is described as in previous work (Lo and Wang, 2006). The SC network generates all-or-none burst activity due to strong recurrent connection there and sends feedback to the Cx, which, due to recruitment of the inhibitory population Cxi, terminates the ramping in the cortical circuit.
The connection efficacies within the BG used in the model are the following for the normal state: gCD–CD = 1.0 nS, gCD–SNr = 3.0 nS, gGPe–GPe = 1.5 nS, gGPe–STN = 0.6 nS, gGPe–SNr = 0.08 nS, gSTN−GPeAMPA = 0.05 nS, gSTN−GPeNMDA = 2.0 nS, gSTN–SNr = 0.06 nS. The efficacy of the connection from the Cx to the CD (gCx–CD) varies between 1 nS and 4.5 nS. The CD to the GPe efficacy gCD–GPe varies between 0 and 8 nS. For the cortical and SC networks, we use the same synaptic connection properties as in past studies (Wang, 2002; Lo and Wang, 2006).
We also investigated how the normal state is transformed into a parkinsonian state characterized by enhanced beta oscillations. Our beta oscillations are generated within the STN–GPe subcircuit. To simulate the parkinsonian state, we update the synaptic efficacies within the STN–GPe circuit to the following values: gGPe–GPe = 0.02 nS, gGPe–STN = 10 nS, gSTN−GPeNMDA = 10 nS (the other connection efficacies within the model are kept the same as the normal state). These updates reflect changes in strengths and numbers of connections as well as lessening of STD of GPe–STN connections resulting from the lower GPe firing rates seen experimentally with dopamine depletion (Stanford and Cooper, 1999; Ogura and Kita, 2000; Cragg et al., 2004; Shen and Johnson, 2005; Kita, 2007; Fan et al., 2012; Wilson, 2013). These and other changes in values were chosen to reproduce key effects qualitatively rather than quantitatively. They amount to a substantial increase of the excitatory–inhibitory interplay between the STN and GPe along with a reduction in the intra-GPe inhibition.
Background inputs, which are incorporated into Equation 4, are described as Poissonian spike trains. All of the BG nuclei in the model receive AMPAergic background inputs, with possible cortical and thalamic origins. GPe neurons also receive GABAergic background inputs representing additional inputs from the striatum. For normal states, the synaptic efficacies and Poisson rates used for the background inputs in our model are as follows: 4.0 nS and 0.8 kHz for the CD; 3.0 nS and 3.2 kHz for GPe AMPA receptors; 2.0 nS and 2.0 kHz for GPe GABAA receptors; 1.6 nS and 4.0 kHz for the STN; 14 nS and 0.8 kHz for the SNr. For the parkinsonion state we use 1.6 nS and 4.0 kHz for GPe AMPA receptors, 22 nS and 2.0 kHz for GPe GABAA receptors, 1.0 nS and 3.2 kHz for the STN, and keep the other values unchanged. The strongly enhanced background CD to GPe efficacy in part represents the removal of the inhibitory effect of dopamine on the striatal neurons projecting to the indirect pathway (Kish et al., 1999; Obeso et al., 2000; Kita, 2007). Furthermore, since differences in ramping CD activity during evidence accumulation between the normal and parkinsonian states have not yet been quantified, we maintain the same ramping pattern across these states and use the enhanced background CD-to-GPe efficacy to represent the overall increase in CD firing rate in the parkinsonian state. If we instead allow the rate of background inputs to increase, then a much smaller increase in CD-to-GPe efficacy yields the same results as we present here (e.g., 5 nS and 10 kHz; data not shown).
Results
Balance of the direct and indirect pathways during evidence accumulation
Previous work suggests that the corticostriatal pathway can effectively adjust the decision threshold during perceptual decision making (Lo and Wang, 2006). The mechanism involved works best in the presence of a strong nonlinearity, i.e., an abrupt suppression of SNr activity immediately before the onset of a saccadic response (presumably coinciding with the threshold crossing), which was observed in physiological recording of monkeys performing saccade tasks (Hikosaka and Wurtz, 1983; Sato and Hikosaka, 2002). With a gradual intensification of CD activity during evidence accumulation in the random-dot task (Ding and Gold, 2010), however, the direct pathway projection would cause the SNr activity to ramp down in a graded way, rather than exhibiting a sudden drop of activity just before a saccade.
By incorporating the indirect pathway as well, we find that SNr activity can remain unchanged over time during most of the CD ramping process (Fig. 2). In our simulations, ramping neuronal activity in the CD leads to a gradual decrease of GPe activity and a corresponding increase of STN activity, which cancels out the direct impact of CD output on SNr activity. This balancing effect weakens, however, when GPe activity becomes sufficiently low and disappears entirely when the GPe activity reaches zero and STN activity correspondingly plateaus. The continued rise in CD output then sharply suppresses SNr activity. In a way similar to that of the earlier model (Lo and Wang, 2006), the abrupt decline of SNr activity releases the downstream suppression of the SC, which combined with direct excitation of decision neurons in the cortex results in a burst of activity in the SC (with peak activity at 200 ms as indicated in Fig. 2F, in this example). The feedback from the SC to the cortex in turn recruits the inhibitory cortical population, Cxi, which terminates the evidence accumulation process there. Thus, the extension of a previous model (Lo and Wang, 2006) to include the indirect pathway allows appropriate evidence accumulation thresholding and decision responses in the presence of more realistic neuronal activity at the striatum. We next consider how this performance depends on the network features that contribute to the balance between the direct and indirect pathway inputs to the SNr.
Behavior of the model Cx–BG–SC circuit when the indirect pathway balances with the direct pathway. A, The two selective populations in the Cx compete with each other and accumulate evidence. B, CD populations innervated by the Cx show ramping activity. C, D, GPe neurons (C) inhibited by CD neurons decrease their activity, while STN neurons (D) increase their activity due to reduced inhibition from the GPe. E, The innervations from the direct and indirect pathways to the SNr neurons balance with each other, such that SNr activity remains unchanged despite the ramping in the CD (blue and green curves). When the activity of GPe neurons decreases to near zero, a continued increase of CD firing leads to further inhibition of the SNr, whose neural activity exhibits a sharp suppression. F, The resulting disinhibition of SC neurons, combined with sufficient cortical excitation, eventually triggers a burst response (indicating a categorical choice) in one of the two SC neural populations (with peak activity at ∼200 ms as indicated by the black bar). Ramping activity in the cortex is terminated by the feedback signal from the SC, defining a threshold for the decision process. Parameter values: gCx–CD = 3 nS, gCD–GPe = 4 nS; coherence level c = 25.6%. The SNr activity in the case when the indirect pathway overbalances the direct pathway, with gCD–GPe = 6 nS, is also shown (E, cyan, left selective; magenta, right selective).
Impact of relative weights of the direct and indirect pathways
We find that when the coherence level of the random-dot sensory inputs is varied, the mean response time increases with the task difficulty (lower coherence), but the decision threshold, defined as the firing rate of the corresponding cortical population at the onset time of SC activity, remains the same (Fig. 3A,B). This finding extends an earlier observation (Lo and Wang, 2006) to the case when gradual evidence accumulation in the CD and the indirect pathway are included. To understand more fully the impact of the indirect pathway, we adjust certain key model parameters and investigate their effects on the decision time and threshold. Adjusting the corticostriatal efficacy, gCx–CD, can effectively modulate the decision threshold, as found previously (Lo and Wang, 2006). Without including the indirect pathway, however, adjusting gCx–CD can only change the threshold value within a narrow range when the ramping activity in the CD during evidence accumulation is taken into account (Fig. 3C, the gCD–GPe = 0 case). We find that enhancing the strength of the indirect pathway by increasing gCD–GPe significantly widens the range of threshold values (Fig. 3C). Defining this range Δth by the difference in threshold values found for gCx–CD between 1.0 and 4.5 nS, we observed a rapid increase in Δth when gCD–GPe is increased from 0 (Fig. 4A). Moreover, Δth saturates when gCD–GPe is larger than a critical value of ∼4 nS (when there is a balance between the direct and indirect pathways; Fig. 2), above which the indirect pathway effectively counterbalances the direct pathway until GPe activity becomes largely suppressed.
Decision threshold and its dependence on the strength of the indirect pathway. A, B, The distribution of decision times and thresholds for two motion-coherence levels. The decision time is modulated by coherence level of the stimulus (top), while the threshold is insensitive to the stimulus strength (bottom). C, The dependence of threshold on the corticostriatal connection strength gCx–CD for different values of gCD–GPe, the strength of the indirect pathway. Error bars indicate SD.
Dependence of range of threshold values on the strength of the indirect pathway. The range of threshold values obtained through gCx–CD modulation (Δth) is enlarged by increasing gCD–GPe. A, When D1-expressing and D2-expressing CD populations have the same firing rates as in Figure 2, the ramp saturates when gCD–GPe exceeds 4 nS, the value for which the indirect pathway balances the direct pathway. B, When D1-expressing and D2-expressing CD populations have different firing rates, the balance value for gCD–GPe increases to 6 nS, since in the present case the D2-expressing CD population selective for motion direction has a lower rate than that of the D1-expressing CD population, and therefore a higher efficacy is needed to achieve balance between the two pathways. Error bars indicate SD.
Striatal projection neurons that mediate direct and indirect pathways express D1 and D2 receptors, respectively (Calabresi et al., 2014). We have assumed that the direct and indirect pathway CD populations have the same firing rates and have lumped them together. The CD neurons projecting through the direct pathway (expressing D1 receptors) and through the indirect pathway (expressing D2 receptors) show dissociable encoding of choice and learning, however (Hikida et al., 2010; Kravitz et al., 2012; Tai et al., 2012). To test the robustness of our decision-making mechanism, we perform additional simulations for the case when left-selective and right-selective D1-expressing and D2-expressing CD neurons have different firing rates, i.e., the CD population selective for the coherent motion direction (say, the left direction) has a higher firing rate than that selective for the right direction for D1-expressing CD neurons, while for D2-expressing CD neurons, the CD population selective for the right has a higher firing rate than that selective for the left direction. This asymmetry of firing rates might result from reward history and, during the preparatory period, will position the subjects to access higher reward opportunities. To simulate our model in this situation, we increased the background efficacy from 4 to 5 nS for D1-expressing CD neurons selective for the left direction and D2-expressing CD neurons selective for the right direction. This introduces an asymmetry in CD firing rates. The resulting network activity resembles that of Figure 2 (data not shown). From Figure 4B, we see that our proposed mechanism about balance between the two pathways still holds. Note that the saturation value for gCD–GPe increases to 6 from 4 nS in our original model (Fig. 4A), since in the present case the D2-expressing CD population selective for left has a lower rate than the rate of the D1-expressing CD population selective for left, and therefore a higher efficacy is needed to achieve balance between the two pathways. This result suggests that the mechanism we have introduced does not depend critically on the assumption that the CD populations projecting to the two pathways have the same spike rates. In the following we will not distinguish D1-expressing and D2-expressing CD neurons since it is currently unknown whether they exhibit different ramping activity during a perceptual decision process.
One possible test of the pathway balancing that we propose would come from localized cooling of particular BG areas, to the point where neuronal responses are slowed but not inactivated (Long and Fee, 2008). Simulations of such cooling in our model, implemented by increasing the neuronal membrane time constant, reveal different effects, depending on the modulated target. Cooling of the GPe has little impact on the balance between pathways, resulting in little change in threshold for particular gCx–CD values and in the overall range of response thresholds (Fig. 5A). In contrast, cooling of the STN weakens the ability of its outputs to counterbalance direct pathway activity. This change results in lowered thresholds over all gCx–CD values as well as a diminished range of response thresholds (Fig. 5B). These alterations and differences between nuclei represent predictions for future experimental investigation.
Impact of slowing down the membrane dynamics of indirect pathway nuclei on threshold modulation by corticostriatal efficacy. A, The dependence on corticostriatal connection strength gCx–CD for different GPe membrane time constants. The time constant τ = C/gL is slowed by increasing the capacitance C with a fixing leak conductance gL. Three time constants are compared, with τ = 20 ms representing the normal time constant of GPe neurons. The threshold is insensitive to the slowing down of GPe membrane dynamics. B, Same as A, but for the STN. τ = 20 ms is the normal time constant of STN neurons. When the STN membrane dynamics are slowed down, the threshold is decreased and its range is reduced.
In our model, the performance and mean decision time (defined as the time lapse from stimulus initiation to threshold crossing) are also influenced by gCD–GPe (Fig. 6). Specifically, when the relative weight of the indirect pathway is enhanced by increasing gCD–GPe from 0 to the balancing value 4 nS, performance improves and response time increases across a wide range of stimulus coherence levels; that is, accuracy is favored over speed. Increasing gCD–GPe even higher, from 4 to 6 nS, yields no change of the performance and mean decision time, consistent with the observation that the decision threshold becomes insensitive to gCD–GPe in this regime of indirect pathway predominance.
Effect of pathway balance on decision making. A, B, Performance (A; proportion of trials for correct response) and mean decision time (B) as a function of stimulus coherence level, with four different gCD–GPe values (gCx–CD = 3 nS). With the strengthening of the input to the indirect pathway, the performance is enhanced, while the decision time is also increased. When the indirect pathway balances or overweights the direct pathway, the performance and mean decision time become saturated (curves with circles and triangles are indistinguishable).
The activity of CD neurons was observed to eventually saturate during evidence accumulation (Ding and Gold, 2010). One way to introduce saturation for CD activity is through the inclusion of STD at corticostriatal synapses. Figure 7A illustrates the mean activation level of AMPAergic synapses with STD receiving a linear ramping presynaptic input. The synaptic activation level Ds, where D and s are shown in the middle panels of Figure 7A, reaches a plateau (Fig. 7A, bottom) depending on the slope of ramping input (Fig. 7A, top). CD activity thus saturates in our circuit model when there is STD at corticostriatal synapses (Fig. 7B). The ramping in the Cx (Fig. 7B, top) leads to a saturation of the activation level of AMPAergic corticostriatal synapses (Fig. 7B, middle), and therefore a saturation of the CD activity (Fig. 7B, bottom). Note that now the input strength to the SNr is weakened due to the CD saturation, and SNr activity remains very high, which would lead to a failure of disinhibition before saccade onset and prolong the decision time. Simply doubling the synaptic efficacies along both the direct and indirect pathways, however, restores the action selection mechanism and the dependence of threshold on gCD–GPe. The saturation of the ensemble-averaged CD activity when aligned with stimulus onset (top) and saccade onset (bottom) are shown in Figure 7C for several stimulus coherence values. The range of threshold values, Δth, as a function of gCD–GPe (Fig. 7D) again shows a two-phase behavior (rapid increase followed by a saturation), resembling that for the original model (Fig. 4A). The critical value of gCD–GPe, above which the indirect pathway balances with the direct pathway, shifts to ∼8 nS now, corresponding to the doubling of efficacies along the direct and indirect pathways. Therefore, taking into account CD saturation does not impair the mechanism that we have introduced for effective threshold modulation.
Decision threshold and its dependence on the strength of the indirect pathway with CD activity saturation. A, Illustration of activation level saturation of an AMPAergic synapse with STD. Top, Linearly ramping presynaptic input rates. Three colored lines correspond to three ramping slopes. Middle, Depression factor D and AMPAergic synaptic gating variable s as a function of time, obtained from Equations 7 and 5 by replacing 
with a population-averaged input rate as shown in the top panel. Bottom, The activation level of an AMPAergic synapse with depression, Ds, shows a saturation depending on the slope of the ramping input. B, Cx and CD activities with the inclusion of STD at corticostriatal synapses in the circuit model. Top, Ramping Cx activity. Middle, Averaged individual AMPAergic synaptic activation level shows a saturation, where Ds = 
with NCx representing the number of neurons in a selective Cx population. Bottom, CD firing rate saturates when the Cx ramps up. C, CD activity for six coherence levels each averaged over 800 trials shows saturation when the activity levels for different trials are aligned with stimulus onset (top) and saccade onset (bottom). D, The range of threshold values, Δth, shows a similar dependence on gCD–GPe as that for original model (compare with Fig. 4A), and saturates when gCD–GPe exceeds 8 nS, the value at which the indirect pathway balances the direct pathway in this case. Error bars indicate SD. Note that the synaptic efficacies along the direct and indirect pathways are doubled to achieve a disinhibition at the SNr before saccade onset despite CD saturation. In B and D, the coherence level is 25.6%. In B and C, gCx–CD = 3 nS and gCD–GPe = 8 nS.
Beta oscillations generated in the STN–GPe subcircuit
Enhanced beta-band oscillations in neuronal population activity in the BG represent one characteristic feature of parkinsonism, which is associated with dopamine depletion in the BG (Boraud et al., 2005; Hammond et al., 2007; Rubin et al., 2012; Calabresi et al., 2014). We investigated how experimentally observed changes induced by dopamine depletion might alter the model behavior, in particular considering the effects of (1) enhanced synaptic connections between STN and GPe neurons (Cragg et al., 2004; Shen and Johnson, 2005; Kita, 2007; Fan et al., 2012; Wilson, 2013), including the STN-to-GPe excitation mediated by NMDA receptors (Bernard and Bolam, 1998; Kita et al., 2004), and (2) enhanced efficacy of the striatal GABAergic projection to the GPe (Kish et al., 1999; Obeso et al., 2000; Kita, 2007).
With appropriate parameter changes reflecting these conditions, the STN–GPe subcircuit of our model gives rise to synchronized firing at a frequency within the beta band (Fig. 8). In this regime, the collective population activity shows an oscillatory pattern (Fig. 8C), with a rhythmic frequency of ∼16 Hz (Fig. 8D). Both STN and GPe neurons exhibit bursts of spikes and also a small proportion of tonic spikes (Fig. 8A,B), emulating the enhanced bursting that can arise in these neurons in parkinsonism (Filion and Tremblay, 1991; Bergman et al., 1994; Wichmann and Soares, 2006). This bursting in our model STN and GPe cells results from their intrinsic PIR property (a depolarizing potential that takes place following prolonged hyperpolarization; Selverston and Moulins, 1985; Llinás, 1988; Getting, 1989; Destexhe and Sejnowski, 2003). In contrast to irregular tonic spiking observed in STN and GPe neurons with control parameter values (Fig. 2, except the rebound bursts of GPe neurons after saccade initiation), enhanced inhibition from the striatum to the GPe and from the GPe to the STN in the model provides the prolonged hyperpolarization necessary for the generation of burst discharges in STN and GPe neurons. Once the T current is available to support PIR, excitatory inputs from the STN can induce bursting in the GPe and excitatory background inputs can help trigger bursting in the STN. While the extent of bursting in our model is exaggerated relative to most experimental conditions, our results nonetheless support the idea that changes in effective connectivity properties associated with parkinsonian conditions can give rise to both beta oscillations, through altered network interactions, and increased bursting, through enhanced recruitment of intrinsic currents.
Beta oscillations in the isolated STN–GPe subcircuit. A, Representative time courses of individual model STN (blue) and GPe (green) neurons. B, Raster plots of model STN (top) and GPe (bottom) neurons (200 neurons of each type). Both STN and GPe neurons exhibit bursting activity, as observed experimentally. C, Population firing rates of STN and GPe neurons. D, Power spectrum of STN and GPe population firing activity, showing a peak around 16 Hz. The prominent beta-range oscillations take place when the GABAA receptor-mediated GPe-to-STN inhibition and NMDA receptor-mediated STN-to-GPe excitation are substantially increased from the control parameter set (see Materials and Methods).
Parameter dependence of beta oscillations and oscillation frequency
Intuitively, the oscillatory periodicity of ∼60 ms arises from the cyclic interplay between NMDA receptor-mediated STN-to-GPe excitation and GABAA receptor-mediated GPe-to-STN inhibition. To better quantify the dependence of the emergence of population oscillations on network properties, we systematically varied several key model parameters. In Figure 9, the effects of varying strengths of connections between the STN and GPe and of background inputs are considered. These variations reveal the parameter regime where the network exhibits beta oscillations (Fig. 9A, top, pale blue region), on which we focus our further analysis. Furthermore, the frequency of beta oscillations is the same in the STN and GPe, reflecting the network-based mechanism underlying this activity. Within the beta oscillation regime, increasing gGPe–STN has little effect on STN firing rate but causes the coefficient of variation (CV) of STN firing to increase. In contrast, increasing gSTN−GPeNMDA strongly decreases the STN firing rate. Together, these two parameters control the GPe firing rate, suggesting that changes in one of these efficacies can counterbalance the other in terms of their influence on this property. Overall, increases in GPe firing do accompany decreases in STN firing rate, because more sustained GPe inhibition provides fewer opportunities for STN bursting.
Dependence of the STN–GPe subcircuit behavior on the interconnection efficacies between the STN and GPe and the background inputs. A, Population oscillation frequency, mean neuronal firing rate, and CV for activity of STN neurons (left column) and GPe neurons (right column) as a function of interconnection strengths gGPe–STN and gSTN−GPeNMDA. Beta-band oscillations are robust over a broad parameter regime. B, Oscillation frequency, firing rate, and CV of STN neurons (left column) and GPe neurons (right column) as a function of the strength of inhibitory background input to GPe neurons (gGPeb) and excitatory background input to STN neurons (gSTNb). In an intermediate range of gGPeb values, CV is large (>0.5; bottom) and the oscillation frequency is in the beta range (∼8–20 Hz; top) for a wide range of gGPeb values. Note that for the yellow regime in the top panels in A, and the red regime in the top panels in B, the CV is very low, and hence we do not consider these as representative of parkinsonian conditions. In determining the population oscillation frequency, we use the autocorrelation function of the population firing rate to estimate the oscillation period T, and the corresponding frequency is given by 1/T.
The GABAergic input from the striatum to the GPe has been suggested to play an important role in determining the oscillation properties in the STN–GPe subcircuit (Terman et al., 2002; Kumar et al., 2011). Our model supports this idea and shows that the background inhibitory input to the GPe, representing input from the striatum and modeled as a background GABAergic synapse with efficacy gGPeb needs to be much stronger than in the normal state for the emergence of beta oscillations (Fig. 2, 2 nS; Fig. 8, 22 nS). Furthermore, we find that the excitatory input to the STN, representing input from the cortex and thalamus and modeled as a background AMPAergic synapse with efficacy gSTNb in our model, can effectively modulate the oscillation frequency within the beta band (high-CV regimes in Fig. 9B, top, ranging from pale blue to yellow). Note that when either gSTNb or gGPeb becomes very large, the oscillation is effectively abolished (indicated by a very low CV in Fig. 9B, bottom). When gGPeb is fixed and gSTNb becomes large (Fig. 9B, top, red regime), the beta oscillation is abolished through a step-like transition from an oscillatory state to an asynchronous state. This transition is induced by the sudden loss of bursting capability in GPe neurons (from bursting dominant to tonic firing dominant) when the excitatory input from the STN becomes too strong, preventing recruitment of the low-threshold T-type calcium current. On the contrary, when gSTNb is fixed (say, 0.8 nS) and gGPeb becomes large (the sector with CV above baseline in Fig. 9B, bottom), the CV decreases gradually as the firing of networks becomes more and more asynchronous, reflecting the stronger random inhibition. The increase of striatal input also decreases the GPe firing rate and correspondingly allows the STN to fire more rapidly (Fig. 9B, middle). The vertical bar in the middle panel of Figure 9B represents a regime where GPe neurons show a mixture of bursts and tonic spikes and STN neurons are silent.
We also investigated the influence of the longest time constants in the network and some other synaptic efficacies on the oscillation properties (Fig. 10). There are three long time constants in our STN–GPe circuit: the NMDA receptor decay time constant, τNMDA, and the recovery time constants of the inactivation variables of the T-type calcium currents responsible for PIR in the STN and GPe, τh+ (STN) and τh+ (GPe), respectively. We see that the oscillation frequencies are insensitive to any individual variation of these long time constants within a wide range, although they have a moderate influence on the firing rates and CVs in the STN and GPe (Fig. 10A–C). A step-like transition occurs when τh+ (GPe) is larger than a critical value (Fig. 10C), induced by a sudden loss of bursting capability of GPe neurons when the recovery of the T-type current in the GPe takes so long that it consistently is interrupted by the excitatory input from the STN, which pulls the GPe neurons out of the bursting regime. The synaptic pathway from the STN to the GPe is mediated through both NMDA and AMPA receptors (Kita et al., 2004). We already showed that when the NMDA component of coupling between the two nuclei is strong, beta oscillations can emerge and persist over a large parameter regime (Fig. 9A,B). The AMPAergic synapse does not influence the oscillation properties when the efficacy is kept at a low level (<1 nS here; Fig. 10D). But if the AMPA strength exceeds a critical value, the beta oscillation is sharply abolished as indicated by a very low CV. This step-like transition occurs when the pulse-like AMPAergic input from the STN is strong enough to interrupt the recovery of the T-type current in the GPe, similarly to the step in Figure 10C. Our model therefore suggests that in the parkinsonian state, the presence of enhanced beta oscillations signifies that the potentiation of the STN-to-GPe projection is predominantly mediated through the NMDAergic synapses.
Parameter dependence of beta oscillations. A–E, Dependence of population oscillation frequency (top), mean neuronal firing rate (middle), and CV (bottom) on the decay constant of the NMDA receptors in the excitatory synaptic projection from the STN to the GPe (A), the slow time constant τh+ of the intrinsic T-type calcium current responsible for PIR in STN neurons (B) and GPe neurons (C), the efficacy of AMPAergic synapses from the STN to the GPe (D), and the strength of recurrent inhibition within the GPe (E). The beta oscillation frequency is insensitive to the exact values of the long time constants (A–C). The step-like reduction of CV in C and D corresponds to a transition from an oscillatory state to an asynchronous state. These transitions are induced by the loss of bursting of GPe neurons due to either an overly long recovery time for PIR (C) or an interruption of the recovery of PIR by an overly strong pulse-like AMPAergic postsynaptic current. An increase of the strength of intra-GPe inhibition significantly reduces CV during beta oscillations (E).
The recurrent connection strength within the GPe, gGPe–GPe, weakly affects the oscillation frequency and firing rate but can more strongly and gradually reduce CV (Fig. 10E). To have a relative large CV in our model, gGPe–GPe needs to be weak (in Figs. 2 and 8, the values of gGPe–GPe used are 1.5 and 0.02 nS, respectively). Recurrent connections within the GPe may be weakened under parkinsonian conditions (Stanford and Cooper, 1999; Ogura and Kita, 2000), but there is conflicting evidence (Miguelez et al., 2012).
Impact of beta oscillations on perceptual decision making
Now we consider the full network model in a parkinsonian state characterized by strong beta oscillations originating in the STN–GPe subcircuit (Fig. 11). Note that when the STN–GPe subcircuit is put back into the full circuit, the GPe receives both background inhibition from the CD and inhibition from the selective populations in the CD and exhibits beta oscillations over a larger parameter regime than in the isolated subcircuit (data not shown). In a simulation of perceptual decision making, neural activity in the Cx and CD modules does not change from the simulated normal state to the parkinsonian one, by the design of the model (compare Figs. 11A,B, 2A,B). The firing rates of the STN and GPe, on the other hand, now oscillate (compare Figs. 11C,D, 2C,D). In the presence of this beta oscillation, as CD activity ramps up over time, the population activity of GPe neurons maintains its oscillatory activity and appears to be insensitive to the CD ramping (Fig. 11C). Oscillations in the STN–GPe subcircuit are transmitted to the SNr (Fig. 11E), and the lack of GPe response to ramping in the CD translates into a loss of the ability of the indirect pathway to balance the direct pathway. Troughs in STN output provide repeated windows of opportunity for the direct pathway to suppress the SNr and allow an SC response (with peak activity at 180 ms, indicated by the black bar in Fig. 11F, in this example).
Decision-making behavior of the full circuit in the simulated parkinsonian state. A, B, Neurons in the Cx (A) and CD (B) show ramping activity, which is similar to the normal case (Fig. 2). C–E, GPe (C) and STN (D) neurons exhibit strong beta-range oscillatory activity, which is transmitted to the SNr (E). F, When the Cx activity is increased to a sufficiently large level and an oscillation results in transient SNr suppression, SC neurons produce a burst of spikes (with peak activity at ∼180 ms as indicated by the black bar) and terminate the ramping activity in the Cx. Parameter values: gCx–CD = 3 nS, gCD–GPe = 4 nS, coherence level c = 25.6%.
As a result, both the decision time and threshold drop substantially in the presence of beta oscillations (Fig. 12). Furthermore, there is a strong reduction in the range of threshold values for the parkinsonian state (Fig. 12C,D), likely because even with increases in gCx–CD, the lack of sensitivity of the beta oscillations to the CD ramping means that qualitatively similar troughs of STN output occur during beta oscillations, again leading to windows of opportunity for disinhibition of SC neurons. Indeed, we find that the decision threshold in simulated parkinsonian states is quite insensitive to beta oscillation frequency, because for all frequencies, repeated epochs arise in which indirect pathway outputs fail to balance the direct pathway (Fig. 12C,D). The speeding up of decision time and loss of threshold variability that we observe are consistent with the accelerated responses observed in parkinsonism using a probabilistic selection task (Frank et al., 2007) and a cognitive task-switching task (Rowe et al., 2008), relative to the control state, and suggest that responses will be less flexible as well, which has been observed previously for Parkinson's disease patients untreated or in the off-medication state in tasks using different cognitive functions (Cools et al., 2001; Monchi et al., 2004).
Reaction time distribution and range of threshold values in the parkinsonian state compared with the balanced normal state. A, B, Reaction-time distribution in the parkinsonian state (green, beta oscillation frequency f = 16 Hz) compared with the normal state (blue). C, D, Threshold variation as a function of the connection efficacy gCx–CD for the normal state (blue) and parkinsonian states with three oscillation frequencies: f = 11, 16, 21 Hz (green). A, C, Coherence level c = 6.4%; B, D, Coherence level c = 25.6%. The decision time is reduced for the parkinsonian state compared with the normal state (A, B), as observed in perceptual decision making in parkinsonism. The range of threshold values is also reduced for the parkinsonian state compared with the normal state (C, D), predicting a loss of flexibility in parkinsonism. The range of threshold values is insensitive to the exact oscillation frequency within the beta band (C, D). In C and D, the background GABAergic conductance to the STN, gSTNb, is tuned to obtain three oscillation frequencies, with gSTNb = 0.7, 1.0, and 2.0 nS, corresponding to f = 11, 16, and 21 Hz, respectively (compare Fig. 9). For the normal and parkinsonian states, gCD–GPe = 4 nS. In A and B, gCx–CD = 3 nS.
Discussion
In this work, we extended a cortico-BG-SC model for threshold detection and modulation (Lo and Wang, 2006) by including the indirect pathway of the BG. Our extended model reproduces the recently observed ramping of CD activity during evidence accumulation in perceptual decision making (Ding and Gold, 2010). The main findings are twofold. First, when the indirect pathway balances the direct pathway, the decision threshold can be modulated by the corticostriatal coupling strength over a wide range, which could be important for flexibly tuning the decision threshold in a speed–accuracy tradeoff. Including STD at corticostriatal synapses in the model reproduces the saturation of CD activity during evidence accumulation (Ding and Gold, 2010) and suggests that the mechanism for threshold modulation still holds when CD activity saturates. Second, we suggest a prominent role for NMDA in generating synchronous beta oscillations that arise in the STN–GPe subnetwork under conditions associated with dopamine depletion, and we predict that in the presence of these oscillations, perceptual decision making becomes more impulsive and less flexible (cf. Cools et al., 2001; Monchi et al., 2004; Frank et al., 2007; Rowe et al., 2008; Milenkova et al., 2011).
A novel feature of our perceptual decision-making model is its emphasis on the GPe contribution to indirect pathway activity. Our model predicts that GPe activity decreases gradually toward zero during evidence accumulation, and the lower bound imposed by the zero activity level provides the nonlinearity needed for a sharp suppression of SNr activity before saccade initiation. This idea could be tested, for example, by unit recordings of GPe neurons in monkeys performing random-dot tasks. Our simulations also predict differential results from slowing time constants of neurons in the GPe and STN, such as by localized cooling that avoids complete inactivation (Long and Fee, 2008). The GPe has previously been suggested to be essential in action selection (Gurney et al., 2001a, b; Frank, 2006; Long and Fee, 2008). In our model, STN activity balances direct inhibition from the CD to the SNr, both delaying reaction time until sufficient evidence has been accumulated and suppressing selection of alternative options. The latter is consistent with a spiking version of an action selection model in which diffuse excitation from the STN to the SNr helps support firing of SNr neurons in nonselected channels; there, however, the STN does not delay direct pathway suppression of the selected SNr channel (Humphries et al., 2006). In the context of decision making, several past works showed an increase or phase shift in STN activity in high-conflict situations (Frank, 2006; Cavanagh et al., 2011; Ratcliff and Frank, 2012; Wiecki and Frank, 2013; Zavala et al., 2013, 2014), possibly due to cortical hyperdirect pathway inputs (cf. Hikosaka and Isoda, 2010). In our model, an increased threshold results from decreased gain of signals from cortex to striatum, and modulations of STN activity, such as by the hyperdirect pathway, could act in addition to this effect, either in a sustained way or to provide rapid decision suppression in response to incongruent information, representing a topic for future work.
The ramping of CD activity following a visual cue may be accompanied by either decreasing or increasing SNr firing rates when monkeys perform a memory-guided saccade task with only one direction rewarded (Sato and Hikosaka, 2002). Our model suggests that these two types of activity patterns might arise when the direct and indirect pathways are balanced and when the indirect pathway overbalances the direct pathway (Fig. 2E, cyan and magenta curves), respectively. The SNr responses could be connected to ramping dopamine signals (Kawagoe et al., 1998, 2004; Hikosaka, 2007), whereas in our model ramping at the CD has a different origin, reflecting a cortical information accumulation process. Ramping striatal activity due to reward accumulation has also been reproduced previously (Ratcliff and Frank, 2012).
Reward-dependent plasticity in the BGs is a form of reinforcement learning that can be incorporated in our model in future research to provide a more realistic implementation of the differential modulation of firing rates of direct and indirect pathway CD neurons by probability of reward, as in Figure 4B. Indeed, an extension of the model with only the direct pathway (Lo and Wang, 2006) that incorporates the dopamine system and reward-dependent plasticity at corticostriatal projections was recently presented (Hsiao and Lo, 2013). Related existing work includes a cortex–BG model implementing optimal decision making for multiple choices (Bogacz and Gurney, 2007) and reinforcement learning (Bogacz and Larsen, 2011) and an extended actor–critic model providing a unified description of reinforcement learning and choice incentive, where the direct and indirect pathway striatal neurons play the roles of two opponent actors and exhibit reward-dependent learning (Collins and Frank, 2014).
The indirect pathway of the BGs is strongly implicated in parkinsonism (Rubin et al., 2012). Parkinsonism is typically associated with delayed motor behaviors; however, experimental results on reaction times of Parkinson's disease patients performing cognitive tasks are highly variable (Robert et al., 2009). These include the observation of decreased reaction times for Parkinson's disease patients off dopamine medication performing a probabilistic selection task (Frank et al., 2007) and a cognitive task-switching task (Rowe et al., 2008), compared with normal subjects, and also the finding that reaction times decrease with STN deep brain stimulation, consistent with the idea that STN modulates decision thresholds (Cavanagh et al., 2011; Green et al., 2013; Frank et al., 2015). In a regime with beta oscillations, our model predicts a reduced decision time for untreated Parkinson's disease subjects performing perceptual decision making. Furthermore, in this regime in our simulations, the range of possible decision thresholds is limited, which predicts that a decreased response flexibility will accompany the reaction-time change, reminiscent of the loss of cognitive flexibility for untreated Parkinson's disease patients in tasks using different cognitive functions (Cools et al., 2001; Monchi et al., 2004). It remains for future work to determine how changes in STN activity would alter decision thresholds in the parkinsonian regime in our model.
Several models about the origins of rhythmic activity, such as beta-band oscillations in parkinsonism, have been proposed, emphasizing the STN–GPe circuit (Terman et al., 2002; Holgado et al., 2010; Kumar et al., 2011; Park et al., 2011; Pavlides et al., 2012), the hyperdirect pathway (Leblois et al., 2006), the striatum (McCarthy et al., 2011), and the cortex (Silberstein et al., 2005; Roopun et al., 2008; Lee et al., 2013). The reciprocal interactions between the STN and GPe provide a natural substrate for rhythm generation (Plenz and Kital, 1999; Bevan et al., 2002b; Terman et al., 2002; Mallet et al., 2008; Tachibana et al., 2011). Explorations of possible mechanisms for beta oscillation rhythmogenesis in the STN–GPe circuit typically involve changes in connection strengths between the STN and GPe or within the GPe, relative to nonoscillatory conditions (Holgado et al., 2010). Our model predicts that potentiated STN-to-GPe connections in parkinsonism are predominantly mediated through NMDA receptors, which can be tested experimentally and used to differentiate this work from previous models. Subunit-specific NMDA antagonists have been found to have antiparkinsonian actions on both motor (Steece-Collier et al., 2000; Löschmann et al., 2004) and cognitive (Hsieh et al., 2012) dysfunctions in animal models of parkinsonism, and the injection of NMDA antagonists directly into the medial pallidum of parkinsonian monkeys has been shown to reverse movement-related symptoms (Graham et al., 1990; Mitchell and Carroll, 1997). Furthermore, consistent with previous work emphasizing the importance of inhibitory input to the GPe (Terman et al., 2002; Kumar et al., 2011), we found that the GABAergic input from the striatum to the GPe, and the excitatory input from the cortex and thalamus to the STN, can modulate oscillation frequencies within the beta band. Importantly, the impact of beta oscillations on perceptual decision making that we demonstrate is independent of the mechanism underlying beta oscillations and would apply even if the beta oscillations were transmitted from a cortical source.
We assume that CD activity in the parkinsonian and normal states is the same, which obviously oversimplifies the neurological changes associated with parkinsonism. The present model can readily incorporate more realistic representations of striatal activity during perceptual decision making for both normal and parkinsonian states when relevant measurements become available. Recently, it was found that the globus pallidus (GP) of rats (homolog of the GPe of primates) is composed of two distinct groups of neurons (Mallet et al., 2008, 2012; Nevado-Holgado et al., 2014), and it would also be of interest to include both types and explore their roles in decision making. Another extension would be the addition of the thalamus to form a complete, spiking cortico-basal ganglia-thalamic loop model. Meanwhile, in its current form, our model suggests novel mechanisms, which should be amenable to future experimental testing, for flexible perceptual decision making in the presence of ramping striatal activity and for how this decision making may be perturbed under the influence of indirect pathway beta oscillations.
Footnotes
This work was supported by a Swartz Foundation Fellowship (W.W.), the Office of Naval Research Grant N00014-13-1-0297 (X.-J.W.), the National Institutes of Health Grants R01 MH062349 (X.-J.W.) and R01 NS070865 (J.E.R.), and the National Science Foundation Division of Mathematical Sciences Award 1312508 (J.E.R.). We thank S. Fusi and C.C. Lo for early work on the computer program used in this study, and the two anonymous reviewers for suggestions in clarifying the paper.
The authors declare no competing financial interests.
- Correspondence should be addressed to Xiao-Jing Wang, Center for Neural Science, New York University, 4 Washington Place, New York, NY 10003. xjwang{at}nyu.edu
