Abstract
Recent physiological studies using behaving monkeys revealed that, in a two-alternative forced-choice visual motion discrimination task, reaction time was correlated with ramping of spike activity of lateral intraparietal cortical neurons. The ramping activity appears to reflect temporal accumulation, on a timescale of hundreds of milliseconds, of sensory evidence before a decision is reached. To elucidate the cellular and circuit basis of such integration times, we developed and investigated a simplified two-variable version of a biophysically realistic cortical network model of decision making. In this model, slow time integration can be achieved robustly if excitatory reverberation is primarily mediated by NMDA receptors; our model with only fast AMPA receptors at recurrent synapses produces decision times that are not comparable with experimental observations. Moreover, we found two distinct modes of network behavior, in which decision computation by winner-take-all competition is instantiated with or without attractor states for working memory. Decision process is closely linked to the local dynamics, in the “decision space” of the system, in the vicinity of an unstable saddle steady state that separates the basins of attraction for the two alternative choices. This picture provides a rigorous and quantitative explanation for the dependence of performance and response time on the degree of task difficulty, and the reason for which reaction times are longer in error trials than in correct trials as observed in the monkey experiment. Our reduced two-variable neural model offers a simple yet biophysically plausible framework for studying perceptual decision making in general.
Introduction
A large body of psychological literature tells us that the time it takes for a choice to be made provides valuable information about decision processes in our mind (Donders, 1868/1969; Posner, 1978; Luce, 1986). Thus, a challenge of considerable interest to neurobiologists is to elucidate the neuronal basis, at the biophysical and circuit levels, of psychophysical reaction times, which are typically many hundreds of milliseconds in nontrivial cognitive tasks. Recently, physiologists using behaving nonhuman primates have begun to reveal firing activities that are correlated with simple decisions (Romo and Salinas, 2001; Schall, 2001; Shadlen and Gold, 2004). In a reaction time paradigm of visual motion discrimination, spike firing of cells in the lateral intraparietal (LIP) cortex of monkeys was found to be correlated with the response time and choice (Roitman and Shadlen, 2002; Huk and Shadlen, 2005). From the onset of a random dot motion stimulus until the time the monkey produces a choice response by a rapid saccadic eye movement, spike activity of LIP neurons selective for a particular saccadic target slowly increases for hundreds of milliseconds. Both this increase in neuronal activity, and the monkey’s behavioral response time, were longer when the percentage of random dots moving coherently (motion strength) was lower. This suggests that LIP neurons could be a candidate system for accumulating uncertain visual information before a perceptual decision is made.
The ramp-to-threshold dynamics is reminiscent of the “diffusion” model (Ratcliff, 1978; Luce, 1986; Smith and Ratcliff, 2004; Palmer et al., 2005), a popular mathematical model used in the study of reaction time tasks. The diffusion model consists of a one-dimensional system that integrates over time the difference between two noisy stimulus inputs. When it reaches one of two thresholds, the choice is made and the decision time is recorded. An important characteristic of the model is that it integrates sensory evidence without any “leakage” (i.e., it is a perfect integrator). The diffusion model fits well to many psychophysical data, is mathematically tractable for analysis, and thus has been a benchmark for other models. Furthermore, it has been shown that the diffusion model can be approximately realized by “connectionist models, ” which may include a leak term; time integration becomes nearly perfect when fine-tuning of parameters cancels out the leakage by network recurrent dynamics (Brown and Holmes, 2001; Usher and McClelland, 2001; Bogacz et al., 2003).
Although the diffusion-type model has also been applied to fit neuronal as well as behavioral data (Shadlen and Newsome, 2001; Mazurek et al., 2003; Ratcliff et al., 2003), its abstract nature does not permit a direct exploration of the cellular and circuit mechanisms that give rise to long integration times in decision processes. In contrast, Wang (2002) investigated a biophysically based cortical microcircuit model for decision making. The model is endowed with slow excitatory reverberation between spiking neurons that produces attractor dynamics, and recurrent feedback inhibition via interneurons underlies winner-take-all behavior. The model replicated most of the psychophysical and physiological results in Shadlen and Newsome (2001) and Roitman and Shadlen (2002).
However, the biologically realistic model consists of thousands of spiking neurons that interact with each other in a highly nonlinear manner. It is difficult to thoroughly analyze the model and understand how it works. For this reason, we have constructed a reduced version of the model in Wang (2002) through a mean-field approach. The simplified model has only two dynamical variables, yet it reproduces much of the behaviors of the original spiking neuron model. The objective of this paper is to present this simplified model and use it to investigate the following biological and conceptual questions. First, a dynamical system with a time constant τ usually can exhibit linear ramping over a timescale limited by τ. Given that the longest biophysical time constant in the model is 100 ms, that of the NMDA receptor (NMDAR)-mediated synaptic current, how does the recurrent dynamics give rise to a much longer integration time τ? Is this slow linear ramping a consequence of a network with slow recurrent excitation? Second, can the model still work when recurrent excitation is solely mediated by the much faster AMPA receptors (AMPARs)? Third, is it necessary that neurons subserving integration during stimulation also show persistent activity during working memory? Finally, what is the relationship, at the mathematical level, between our neuronal model and the diffusion model? Answers to these issues will help to elucidate the cortical circuit mechanisms of perceptual decision process.
Materials and Methods
Model reduction
Simplified mean-field approach
Details of the original network of spiking neurons used for making binary decisions (Wang, 2002) can be found in supplementary information A (available at www.jneurosci.org as supplemental material). Here, we used a “mean-field” approach to reduce the model. This approach has been used to analytically study spiking neuronal network models comprised of integrate-and-fire types (Amit and Brunel, 1997a, b; Brunel and Wang, 2001; Renart et al., 2003). Briefly, the net input to a neuron in a large homogeneous population is treated as a Gaussian random process. Then, it can be shown that the mean activity of a (homogeneous) population can be represented by a single unit (Fig. 1). The population firing rates depend on the input currents, which in turn depend on the firing rates. Thus, the population firing rate (or the firing rate of a representative neuron) must be determined self-consistently. These calculations are computationally extensive, taking into account realistic synaptic dynamics (Wang, 1999; Brunel and Wang, 2001) and higher order corrections (Fourcaud and Brunel, 2002). In this work, we propose to use a more simplified approach. First, the driving force of the synaptic currents are assumed to be constant as in Brunel (2000). Second, we assume that the variance of the membrane potential of the cell is mainly contributed by the external input to each neuron, while the contributions from the recurrent connections are averaged out because of the all-to-all connectivity (Renart et al., 2003) and by the averaging effect of the long time constant of NMDA receptors. The SD of fluctuations σ does not vary significantly and thus are fixed as constant.
Reduction of a biophysical neuronal decision-making model. The original model (top) is endowed with strong recurrent excitation between neurons with similar stimulus selectivity, and effective inhibition between them via shared inhibition. NS and I denote the nonselective excitatory (black) and inhibitory (green) pools of cells, respectively. Arrows, Excitatory connections; circles, inhibitory connections. I1 and I2 are inputs from external stimulus to selective neural populations 1 (blue) and 2 (red). Brown arrows, Background noisy inputs. w+ denotes enhanced excitatory connections within each selective neural pool. The numbers on the right displays the total number of dynamical equations involved in the model. First step, Mean-field approach reduces 2000 spiking neurons into four neural units (with a total of 11 dynamical variables). Second step, Simplify the linear input–output relation (F–I curves) of the cells: (1) fit the input–output relation (F–I curve) of the spiking neuronal model with a simple function (Abbott and Chance, 2005); (2) linearize F–I curve for I cells; and (3) assume constant activity of NS cells. The final step involves the assumption that all fast variables of the system reach steady states earlier than that of NMDAR. The final reduced two-variable model (bottom) consists of two neural units, endowed with self-excitation and effective mutual inhibition.
Moreover, the firing rate r of a leaky integrate-and-fire (LIF) neuron receiving a noisy input current is given by the first-passage time formula (Ricciardi, 1977; Amit and Tsodyks, 1991; Renart et al., 2003): 
(1) where φ is a function of the total synaptic input current Isyn. τm is the membrane time constant. Vth is the spiking threshold for the membrane voltage, Vreset is the reset voltage, τref is the refractory period, σV is the membrane potential SD, and Vss = VL + Isyn/gL. Instead of using this formula, we adopted a simplified input–output function from Abbott and Chance (2005) as follows: 
(2)
In this equation, the subscripts E and I are labels for a pyramidal neuron and an interneuron, respectively. Isyn is the total synaptic input to a single cell, and cE, I is the gain factor. gE, I is a noise factor that determines the shape of the “curvature” of φ. If gE, I is large, φ would act like a linear-threshold function with IE, I/c as the threshold current (supplemental Fig. 1, supplementary information B, available at www.jneurosci.org as supplemental material). These parameters were obtained by fitting the model to the first-passage time formula Equation 1 of a single-cell LIF model driven by AMPA receptor-mediated external Gaussian noise (corresponding to a Poisson input at 2.4 kHz) (supplemental Fig. 1, supplementary information B, available at www.jneurosci.org as supplemental material). The resulting parameter values are, for pyramidal cells, IE = 125 Hz, gE = 0.16 s, and cE = 310(VnC)−1; and for interneurons, II = 177 Hz, gI = 0.087 s, and cI = 615(VnC)−1. The fit by Equation 2 is particularly accurate when the cells receive large noise and achieve moderate firing rates.
With these simplifications, we can reduce the spiking neural network to a system with 11 variables, describing the mean firing rates and the output synaptic gating variables of four different neural populations (two selective and one nonselective excitatory cell populations, and one inhibitory cell population). The mean-field theory is a framework for calculating steady states but does not provide a systematic recipe for describing temporal dynamics. We assume that the population firing rate of each population can be described by the Wilson–Cowan type of equations (Wilson and Cowan, 1972, 1973) with a fast time constant (τr = 2 ms) in noisy networks (van Vreeswijk and Sompolinsky, 1998; Brunel et al., 2001; Fourcaud and Brunel, 2002; Renart et al., 2003). The 11 dynamical equations are as follows:
(3)
(4)
(5)
(6)
(7)
where i = 1, 2, 3 denotes the two selective, and one nonselective excitatory populations, and I is the inhibitory population. ri(t) is the instantaneous mean firing rate of the presynaptic excitatory population i, and rI(t) is the mean firing rate of the inhibitory population. S and its associated τ are the average synaptic gating variable and its corresponding decay time constant, respectively, with their receptor type denoted by their subscripts. F(ψi) = ψi/(τNMDA(1 − ψi)), and ψi is the steady state of Si.
The dynamics of the NMDA gating variable is characterized by a fast rise followed by a slow decay (Wang, 1999). Given that the presynaptic inputs at the recurrent synapses are described by a train of delta-like spikes, and assuming the interspike intervals to be nearly Poisson, the average gating variable can be fitted by a simple function (supplemental Fig. 2, supplementary information C, available at www.jneurosci.org as supplemental material): 
(8) where γ = 0.641 and r is the presynaptic firing rate. It can be easily shown by simple algebra that F(ψ(r)) = γr [but see Brunel and Wang (2001) for a more rigorous treatment].
Time course with two different motion strengths. Motion coherence of 0% (black traces) and 51.2% (light gray traces) each with 10 sample trials. Firing rates that ramp upward (bold traces) are for saccades made toward the RF of the neuron, whereas downward (dashed traces) are for saccades away from RF. Ramping is steeper for higher coherence level. The prescribed threshold is fixed at 15 Hz. Once the firing rate crosses the threshold, a decision is made, and the decision time is the time it takes from stimulus onset (0 ms) until the threshold is crossed. The reaction time is defined as the decision time plus a nondecision latency of 100 ms. The bold horizontal line at the top of the figure denotes the duration, at zero coherence, where the firing rates toward and away from RF are indistinguishable.
Phase-plane reduction
The model can be further reduced to a two-variable system. This is done through the following three approximations.
(1) Constant activity of nonselective cells.
Under a wide range of conditions, the firing rate of the nonselective population changes only by a modest amount. Thus, we assume that the nonselective populations fire at a constant mean rate of 2 Hz. This reduces the system to three populations (Fig. 1). A consequence of this approximation is that the overall excitatory drive is no longer normalized and that the spontaneous state of the nonselective population does not equal that of the selective ones (Amit and Brunel, 1997b). Nevertheless, the difference in spontaneous firing rates between them is small (1 Hz). Another consequence of this assumption is that the extra inhibition on the selective populations contributed by the slightly elevated activity of the nonselective population (through the interneurons) would be neglected.
(2) Linearization of the input–output relation of the interneuron.
In principle, the inhibitory population rate depends on itself (via inhibitory–inhibitory coupling) as well as the excitatory firing rates (via excitatory–inhibitory coupling) and hence is not given explicitly. This complication, however, can be eliminated by a linear approximation of the input–output transfer function of the inhibitory cell. Unlike the excitatory cells, in the network in Wang (2002), the interneurons have a higher spontaneous firing rate of ∼8 Hz. The mean firing rate of the inhibitory population typically falls between the range of 8–15 Hz. Within this range, the single-cell input–output relation is almost linear (supplemental Fig. 1B, supplementary information B, available at www.jneurosci.org as supplemental material), fitted by the following: 
(9) where g2 = 2 and r0 = 11.5 Hz.
With the linearity in the response of the inhibitory population, the self-inhibitory term can now be easily absorbed into the factors g2 and r0. In particular, if we define JII as the self-inhibitory synaptic coupling, then self-inhibition is expressed as a term −JIIφ in Isyn. Grouping the two φ dependent terms together, it is clear that self-inhibition term effectively lowers the mean firing rate of interneurons, rI, by a factor of 1 + (cI/g2)JII. This helps to simplify the calculations by removing self-consistency calculations for the inhibitory population (Fig. 1).
(3) Slow dynamics of NMDA gating variable.
Our model involves membrane time constants of neurons and those of synaptic gating variables. For the LIF neuron model driven by filtered noisy inputs, it has been shown that the firing response to a stimulus is instantaneous (Brunel et al., 2001; Fourcaud and Brunel, 2002). Hence, the membrane time constant of the single cell can be neglected. Furthermore, among the synaptic transmissions mediated by AMPA, NMDA, and GABAA receptors, the synaptic gating variable of NMDA receptors has the longest decay time constant (100 ms). Therefore, we can assume that all other variables achieve their steady states much faster than the NMDA gating variable SNMDA, which dominates the time evolution of the system. Specifically, we define the two dynamical equations of the system as follows (Wang and Rinzel, 1992; Ermentrout, 1994; Renart et al., 2003): 
where i = 1, 2 labels the two excitatory populations (selective for left or right motion directions). The AMPA and GABAA gating variables reach their steady states much faster than that of NMDA receptors, which means that the average gating variables of AMPA and GABAA receptors become proportional to the average firing rates of presynaptic cells (Brunel and Wang, 2001):
Dynamical equations
In summary, we have reduced a network model with two thousand spiking neurons into a two-variable system (Fig. 1) described by the dynamical equations as follows: 
(10) 
(11) where the two excitatory neural populations (selective for rightward and leftward motion directions) are labeled by 1 and 2, and, for the sake of convenience, we denote S for SNMDA and τS for τNMDA. The firing rates r1 and r2 are given by Equation 2: 
(12) 
(13) 
(14) 
(15) where Ii represents the visual motion stimulus to the population i and depends on the motion strength (see Results). Inoise, i is a noise term, and I0 is the mean effective external input common to both populations. Because of the background input into nonselective excitatory cells and interneurons, I0 includes not only direct background input to a selective population but also indirect background inputs from these nonselective cells. Isyn, i is the total synaptic current from both recurrent connections and inputs fed from outside the local network. The coefficients JN, ij and JA, ij are effective coupling constants from neuron j to i mediated by NMDAR and AMPAR, respectively. The negative sign in front of JN, 12, JN, 21, JA, 12, and JA, 21 indicates that the overall effective connectivity between the two selective populations is inhibitory. This is because the inhibitory cell population (I) receives inputs from both excitatory neural populations (E); hence its output (proportional to input) is of the form JN, I→EJE→I(S1 + S2) and JA, I→EJE→I(r1 + r2). Thus, for example, in Isyn, 1, the S2-dependent term is of the form (JN, E→E − JN, I→EJE→I)S2 ≡ −JN, 12S2 with JN, 12 = JN, I→EJE→I − JN, E→E· JN, 21, JA, 12, and JA, 21 are defined similarly.
Note that, in solving Equations 10–15, a practical difficulty resides in the fact that the firing rates are not given explicitly: ri = φ(Isyn, i), where Isyn, i depends on both r1 and r2. This problem was solved by finding an effective single-cell input–output relation H (for details, see supplementary information D, available at www.jneurosci.org as supplemental material). Thus, Equations 12–15 are rewritten as follows: 
(16) 
(17) 
18 
(19)
Combining Equations 16–19 with Equations 10 and 11, we have finally a self-contained two-variable system for S1 and S2: 
(20) 
(21)
Parameter values
In our model, the decision (1 = right; 2 = left) is made when one of the two competing neural populations reaches a fixed threshold, for example, θ = 15 Hz. The decision time is defined as the time it takes for the activity of the “winning” population to travel from its initial (spontaneous) state to the decision threshold θ. The sum of the decision time and a fixed nondecision time constant (reflecting a combination of sensory input latency and motor response), which we chose to be 100 ms (but see Mazurek et al., 2003), yields the reaction time. Unless otherwise mentioned, the standard set of parameters for the two-variable model is as follows: JN, 11 = 0.1561 nA = JN, 22, JN, 12 = 0.0264 nA = JN, 21, JA, 11 = 9.9026 × 10−4 nC = JA, 22, JA, 12 = 6.5177 × 10−5 nA · Hz−1 = JA, 21 and I0 = 0.2346 nA. These values are deduced from the parameters of the original spiking neural network model and are slightly adjusted so that the model reproduces the reaction times observed in the monkey’s experiment (Roitman and Shadlen, 2002). (Note that this set of parameters may not be optimal in fitting the experimental data.)
Simulations
The mean-field approach does not include time-varying noise that plays a critical role in the spiking neural network. To amend this, we added a noise term Inoise implemented as a white noise filtered by a short (AMPA synaptic) time constant. This is thus described by an Ornstein–Uhlenbeck process (Uhlenbeck and Ornstein, 1930) (for example, see Destexhe et al., 2001): 
where σnoise2 is the variance of the noise, and η is a Gaussian white noise with zero mean and unit variance. Unless specified, σnoise is fixed at 0.007 nA.
Simulation codes were written in Matlab, run on a Linux workstation. Because the reduced rate model does not require a high temporal resolution, Euler’s method with an integration time step of 0.1 ms was used for numerical integration of the dynamical equations. Simulation results were checked and confirmed using smaller time steps (down to 0.01 ms). The instantaneous population firing rates were calculated by averaging over a time window of 50 ms, slided with a time step of 5 ms. For computation of the psychometric and chronometric functions, each data point was obtained from a block of 2000 trials.
Phase-plane and bifurcation analyses of the reduced model were performed using XPPAUT (Ermentrout, 1990).
Results
Comparison between model and experiment
The inputs to our LIP model neurons should mimic the output of upstream neurons in middle temporal (MT) area that encode the visual motion stimulus. Following Britten et al. (1993) and Wang (2002), the firing rate of an MT cell is a linear function of the motion strength (percentage coherence) c′, increasing or decreasing with c′ depending on whether the motion stimulus is in the preferred or nonpreferred (null) direction of the cell. The (absolute) stimulus strength, μ0, is the input received when there is no bias (c′ = 0%). Note that we do not include noise in the stimulus input because it was found that trial-by-trial stimulus variability is not the primary source of stochasticity in decision choices (Britten et al., 1996).
Specifically, if the stimulus is biased, favoring population 1 with a nonzero coherence level c′, then the synaptic currents attributable to the stimulus alone are as follows: 
where JA, ext = 0.2243 × 10−3 nA · Hz−1 is the average synaptic coupling with AMPARs.
In Figure 2, we show sample time course for two different coherence levels c′. We can see that as c′ increases, the ramping activity of the neural population, whose response field (RF) is the saccadic target, becomes steeper. Therefore, the decision time is shorter for higher c′. This is an expected result because the higher the overall external inputs, the steeper will be the ramping activities, and the faster will be the accumulation of sensory evidence before a decision is made. Note that during motion viewing of lower coherence levels, there is an initial time epoch lasting for hundreds of milliseconds (denoted by a black horizontal bar) when the two traces of activity are indistinguishable before they eventually split apart. This biphasic phenomenon will be explained in later analysis.
The two-variable model replicates fairly well the psychometric function (behavioral performance) and chronometric function (reaction time of correct and error trials) of both the original large-scale spiking neuronal network model (Wang, 2002) and the monkey experiment (Roitman and Shadlen, 2002) (Fig. 3). The psychometric and neurometric functions in Figure 3 are fitted with a Weibull function (Quick, 1974): 
where p is the probability of a correct choice, α is the discrimination threshold at which the performance is 82% correct, and β describes the slope of the psychometric function. With the set of parameters that we used, the reduced model has a threshold α of 7.2% and a slope β of 1.25. These values are comparable with the experimental values of 7.4% and 1.3, and that of the full spiking network model of 8.4% and 1.6.
Performance and reaction time of models and the experiment of Roitman and Shadlen (2002). First column, Psychometric data from experiment and the models (data are fit with a Weibull function). Second column, Reaction time from experiment and the models. Open circles joined by dashed lines, Mean reaction of error trials; filled circles, correct trials. σnoise = 0.008 nA. Experimental data are adapted from Mazurek et al. (2003).
Decision space analysis of time integration and categorical choice
Random decision with unbiased external stimulus
To investigate how the network model responds to, and integrates over time, a stimulus, we performed a phase-plane analysis of the model system. This is done by first setting the dynamical equations dS1/dt = G1(S1, S2) = 0 and dS2/dt = G2(S1, S2) = 0, and then plotting these two lines in the (S1, S2) phase (decision) space. These two lines are called nullclines, and their intersections are steady states of the system. Furthermore, the stability of the steady states is determined by how the nullclines intersect with each other (Strogatz, 2001). Because the average Si of population i is a monotonic, increasing function of the population firing rate ri and their steady states are related by the simple monotonic Equation 8 (see Materials and Methods), we expect the nullclines in the (r1, r2) space to be qualitatively similar.
In the absence of a stimulus, the two nullclines intersect with each other five times, producing five steady states, of which three are stable (attractors) and two are unstable (Fig. 4A). In the absence of stimulation, the network is unbiased and lies at the lower-left symmetrical attractor state, corresponding to the spontaneous state where both populations fire at a low rate. The two (upper left and lower right) asymmetrical attractors correspond to mnemonic persistent states in which one of the neural populations exhibit self-sustained elevated spike activity. Thus, our model can subserve working memory: a transient input can bring the system from the resting state to one of the two stimulus-selective persistent activity states, which can be internally maintained across a delay period.
Random choice with stimulus at zero coherence. A, Phase-plane without stimulus. Black circles, Stable steady states (attractors); gray circles, saddle-type unstable steady states. The green and orange lines are the nullclines for the synaptic dynamical variables S1 and S2. Using Equation 8, a threshold at 15 Hz would correspond to S = 0.49 in phase space. B, With an unbiased stimulus of 30 Hz, the two unstable steady states together with the low stable steady state disappear, and a new symmetric unstable steady state is formed. The black line with arrows toward (away) from the saddle point is the stable (unstable) manifold of this saddle point. The stable manifold is exactly the boundary between the two basins of attraction of the two choice attractors (when there is no noise). Superimposed are two typical single-trial trajectories (blue and red lines) of the state of the system from simulations. Color labeling is the same as in Figure 1. C, Schematic diagram of a generic saddle-like steady state (gray circle) and the local flows (arrows) around it. The lines directly toward (magenta) and away (brown) the steady state (gray) are its stable and unstable eigenvectors, respectively, with an exponential temporal dynamics determined by τstable (brown) and τunstable (magenta). D, A diagram of how a one-dimensional decision “landscape” changes with stimulus inputs in a single trial, illustrating decision computation and working memory by the same network. See Results for detailed description.
When a stimulus (e.g., with μ0 = 30 Hz and c′ = 0%) is applied, the phase space of the model is reconfigured. The spontaneous state vanishes. At the same time, a saddle-type unstable steady state is created that separates the two asymmetrical attractors (Fig. 4B). In Figure 4B, two lines emanate from this saddle point (Fig. 4C). One of them is called the stable manifold. The system starting at any point on this curve eventually converges to the saddle point. The stable manifold forms a boundary that separates the two basins of attraction: if the (noiseless) system starts within a basin of attraction (left or right from the stable manifold), it will be attracted toward the associated asymmetric attractor. The other, unstable, manifold extends from the saddle point to the attractors. The system starting at a point on this manifold is repelled to one of the two competing attractors. Although the phase space is symmetrical, the addition of noise can perturb the system to move away from the diagonal line, and toward one of the two competing attractors, so that a categorical choice is made by the model. Therefore, this picture offers a mathematical description of a two-alternative forced-choice computation, even when at zero coherence the average sensory input is the same for each of the two neural populations. More specifically, as illustrated by the simulation results from two individual trials (Fig. 4B, red and blue lines), the temporal dynamics of decision making consists of two steps: the system initially wanders along the diagonal line (when the two population rates are about the same), before it converges to one of the two asymmetrical attractor states (when one of the populations increases, while the other decreases) corresponding to a categorical choice. Interestingly, there is evidence that LIP neurons recorded from behaving monkey during the visual motion discrimination experiment also exhibit such biphasic time courses (Roitman and Shadlen, 2002; Huk and Shadlen, 2005). This characteristic will also be important when we compare our model with the diffusion model.
In a delayed response version of the task (Shadlen and Newsome, 2001), the motion stimulus is shown for a fixed duration, and the monkey is required to withhold the response across a delay period of a few seconds when the choice must be maintained actively in working memory. Our model is able to perform the delayed response task, because after the stimulus offset the phase-space configuration of the system reverts to that of Figure 4A, in which the choice (1 or 2) is stored in one of the two asymmetric attractor states, and that information can be retrieved at the end of the delay to produce a motor response (leftward or rightward saccade). Another transient input can bring the system back to the resting state, erasing the memory trace. We shall discuss later in more detail the relationship between decision making and working memory.
This phase space depiction is to be contrasted with a schematic view of the decision-making dynamics in terms of a one-dimensional “energy landscape” (Fig. 4D). In it, a hypothetical “energy function” is plotted as function of a single dynamical variable of the system (in our case, S1 − S2 or r1 − r2). The energy always decreases as the system evolves in time, and each local minimum of the energy function is an attractor (Amit, 1992). The local maxima are the unstable saddle points that separate the basins of attractions. The instantaneous state of the system is indicated by a ball. Red, blue, and black portions of the landscape denote the basins of attraction of the two competing attractors and the spontaneous state, respectively. Black arrows denote the most likely direction of motion of the ball. When an unbiased stimulus is presented, the spontaneous state disappears and a saddle point emerges.
This instability, coupled with noise, forces the network to approach one of the attractors. After the stimulus offset, the system can store the choice in short-term memory by virtue of the asymmetrical attractor states (see details below). Although instructive, Figure 4D is purely illustrative. In particular, it assumes that the dynamics is one-dimensional, although we have seen that the dynamics of the model should be understood in the two-dimensional phase space. Below, we will address the question of whether the system can be reduced to a one-dimensional dynamics under certain conditions.
Decision with nonzero coherence: biased basins of attraction
We have seen in Figures 2 and 3 that by increasing c′, we can have a steeper ramping activity and a corresponding shorter decision time. In Figure 5, A and B, we show how the phase space changes when a weak motion stimulus (c′ = 6.4%) is presented. The phase space is no longer symmetrical: the attractor state 1 (correct choice) has a larger basin of attraction than attractor 2. This is because the neural population 1 receives a stronger external input, raising its nullcline (green) higher than that of population 2. At the onset of a biased external input, the spontaneous state of the system already lies within the basin of attraction of the favored choice attractor 1 (for a clearer example at a higher coherence, see Fig. 5C). This leads to a higher percentage of correct choices above the chance level. Shown in Figure 5B are also the dynamical trajectories of the system in one correct trial and one error trial. In an error trial, because the system is initially in the basin of attractor 1, it has to travel across the stable manifold of the saddle to reach the basin of attractor 2. In doing so, the system has a tendency to stay close to the stable manifold and move toward the saddle point, then diverging from it along the unstable manifold (Fig. 5B, trajectory in red). Recall that a saddle point is a steady state, and the dynamics near it is very slow. Specifically, at a distance δ away from the saddle point, the time for the system to stay around that area is ∼log(1/δ) (Hubbard and West, 1995), which goes to infinity as δ approaches zero. (This will become clearer when we discuss the local dynamics near a saddle point.) Thus, the time it takes for the system to reach one of the choice attractors is on average longer in error trials than in correct trials (Fig. 3). This provides a rigorous mathematical explanation as to why in our model the decision times in error trials are generally longer than in correct trials, a salient behavioral observation in the monkey experiments (Roitman and Shadlen, 2002; Mazurek et al., 2003) that is reproduced by our model.
Basins of attraction with stimulus at nonzero coherence (c′ > 0%). A, Phase-plane without stimulus as in Figure 4A. B, The stable manifold is tilted away from the spontaneous state and toward the less favored attractor when c′ is nonzero (6.4%). As a result, at the onset of stimulus, the system starts in a resting state that has a higher chance of falling in the basin of attraction of the favored attractor state. The blue and red lines are typical single-trial trajectories for correct and error choices, respectively. C, Stronger bias between the basins of the two competing attractor states with a larger c′ (=51.2%). D, When c′ is sufficiently large, the saddle steady state annihilates with the less favored attractor, leaving only one choice attractor. c′ = 100%.
As c′ increases, the basin of attraction of the favored attractor increases at the expense of its competitor (Fig. 5B, C). This results in a higher probability of making correct choices. The correct decisions are made more quickly at higher c′ because the spontaneous state (the starting point of the system at the onset of stimulation) lies deeper into the favored attractor’s basin, and the saddle point is remote from the trajectory of the system in correct trials. With a sufficiently large c′ (Fig. 5C, c′ = 51.2%), the distance from the stable manifold to the spontaneous state is so large that fluctuations attributable to noise are insufficient to bring any trajectory across the stable manifold to the less favored attractor. Hence, the system only makes correct choices. Varying the noise level will affect how large c′ must be for the performance to be 100%. However, there is a critical c′ level (∼70%), above which the less favored attractor annihilates with the saddle point, so that only the favored attractor exists and the system always make the correct choice regardless of the noise level (Fig. 5D).
Integration time and local dynamics near the saddle point
A central issue of the present study is to understand how a decision is made at very low or zero motion strength, and what determines the integration time much longer than the biophysical time constants inherent in the model. We have seen in Figures 4 and 5 that, at the onset of a stimulus, the system is in the resting state that is either on the stable manifold (if c′ = 0%) or very close to it (if c′ is small). Therefore, the system evolves toward the saddle point along the stable manifold, before diverging from it (because of noise in the case of c′ = 0%) along the unstable manifold and eventually reaching one of the two choice attractor states. The dynamics in the vicinity of the saddle point is slow, and thus is expected to contribute greatly to the integration time of the system.
Near the saddle point the dynamics of the system can be linearized so that the difference ΔS = S − Ssaddle, between S = (S1, S2) and the saddle Ssaddle can be written as follows: 
(22) where v1 and v2 are the “eigenvectors” of the saddle, which are the same as the tangents of the two manifolds at the saddle point (Strogatz, 2001) (Fig. 4D); −1/τstable and 1/τunstable are the stable and unstable eigenvalues of the saddle. The coefficients a1 and a2 are given by the initial condition of the system. For example, suppose the system starts on the stable manifold toward the saddle point, so that a2 = 0. The distance between the initial state of the system to this steady state after a time interval t decreases as ∼exp(−t/τstable). In contrast, if the system moves along the unstable manifold, a1 = 0, it is repelled away from the saddle with the distance increasing as ∼exp(t/τunstable). In general, a trajectory in a single trial would be a linear combination of these two exponentials (Eq. 22). Note that, if the system is at a distance δ from the saddle, then the time for ΔS to become large (of order 1) is given by 1 ≃ δexp(tstay/τunstable), which yields tstay ≃ τunstablelog(1/δ).
As μ0 is decreased, neurons take longer time to accumulate information from a weaker stimulus, and the decision time increases monotonically (Fig. 6A, B). We found that τunstable of the saddle point increases in parallel with the decision time (Fig. 6C). Note that τunstable diverges to infinity as μ0 approaches a bifurcation point at μ0* ≃ 7 Hz (see Figs. 10 and 6C), but the decision time remains finite (Fig. 6A). Moreover, at small μ0, τunstable is significantly larger than τstable, in which case one can qualitatively reduce the model to a one-dimensional diffusion like system (supplementary information E, available at www.jneurosci.org as supplemental material). Quantitatively, however, this is not a good approximation. To describe the system only in terms of the difference S1 − S2 (along the unstable manifold of the saddle), one must disregard the initial phase of the decision dynamics along the stable manifold, which, with a time constant of ∼200–300 ms, contributes greatly to the decision time and should not be neglected.
Decision time and local dynamics in the vicinity of a saddle point. Zero coherence level from A to C. A, Longer reaction time for smaller stimulus strength μ0. Error bars indicate SD. B, Typical time courses: ramping is faster for larger stimulus strength, μ0. C, Time constants of saddle-like unstable steady state with different μ0. For μ0 > 17 Hz, τstable is larger than τunstable, whereas the opposite is true for μ0 < 17 Hz. D, Time constants of the unstable saddle as function of coherence level c′ (μ0 fixed at 30 Hz). The unstable time constant is essentially constant up to c′ ∼ 70%. The sudden increase in τunstable happens just before the bifurcation point at which the saddle coalesces with the less favored attractor and disappears (see Fig. 5).
In contrast, for a reasonably large μ0 (e.g., 30 Hz), τstable is longer than τunstable (Fig. 6C); therefore, the system definitely cannot be described by a one-dimensional dynamics along the unstable manifold. In fact, as we have seen in Figures 4B and 5B, noise often perturbs the system away from the stable manifold, before the trajectory gets close to the saddle point, so that the system never approaches the unstable manifold of the saddle and τunstable can no longer be critically important. Moreover, because τstable and τunstable are relevant only if the system passes through the vicinity of the saddle point, they have an increasingly weak effect on the decision time with a larger μ0.
This is also true with an increased coherence level c′ (for a fixed μ0 = 30 Hz), in which case the system rarely gets a chance to be close to the saddle point (Fig. 5C). We found that, whereas the decision time of the system monotonically decreases with increasing c′ (Fig. 3), the two characteristic time constants of the saddle point do not vary significantly with c′ (Fig. 6D), except for c′ > 70% where the divergence of τunstable is associated with the disappearance of the saddle point (compare Fig. 5D).
Slow decision time with NMDA receptors
An important finding from the above analysis is that integration time of many hundreds of milliseconds are realized robustly without the need of fine-tuning parameters (such as μ0) (Fig. 6A). Previous simulation results (Wang, 2002) indicated that this desirable feature critically depends on the NMDARs at recurrent excitatory synapses. To assess whether NMDARs are required, it is necessary to incorporate the AMPARs and consider the time integration of the system as the AMPA/NMDA ratio is varied. To this end, we could no longer use the two-variable model, which does not take into account the fast AMPA dynamics. Instead, we investigated the full, 11-variable, dynamical system (Eqs. 3–7).
We calculated the AMPA:NMDA ratio from their relative contributions to the unitary synaptic current at recurrent connections, defined by the charge (time integral of current) elicited at −65 mV by a single presynaptic spike (Compte et al., 2000). With our standard parameter set, the AMPA:NMDA ratio is 15:85. To vary this ratio, the maximum synaptic conductances, gAMPA and gNMDA, at all recurrent excitatory connections were changed at the same time, with the total (summated) unitary charge conserved. (Note that, because the NMDA conductance is voltage dependent and summates over time, the effective AMPA:NMDA ratio varies with the postsynaptic firing rate and thus cannot be fixed by our method.)
We found that, with an increased AMPA:NMDA ratio at the recurrent synapses, the activity becomes much faster (Fig. 7A). Steeper ramping activity results in a shorter decision time (Fig. 7B). At the same time, the decision performance significantly deteriorates (Fig. 7C). When the psychometric function was fitted by a Weibull function, we found that the threshold (α) increases with an increased AMPA:NMDA ratio. Their values are 4.79, 6.11, and 9.37% with 0, 25, and 50% AMPA at the recurrent connections. In contrast, the slope (β) is 1, 1.38, and 1.32 respectively. The shorter reaction time holds true not only as function of the coherence level (Fig. 7B) but also as function of the stimulus amplitude (Fig. 7D).
Dependence of decision-making behavior on the AMPA:NMDA ratio at recurrent synapses. A, Typical time courses: faster ramping neural activity at larger AMPA:NMDA ratios. Top black (gray) horizontal bar denotes the duration where the firing rates are not distinguishable [i.e., the trajectory lies along the stable manifold of the saddle point, when AMPA:NMDA is 35:65 (0:100)]. B, Reaction time is shorter with a higher AMPA:NMDA ratio. C, The performance, however, becomes less accurate. Accuracy data are fitted by a Weibull function. D, For c′ = 0%, a higher AMPA:NMDA ratio decreases the reaction time for the entire range of stimulus strengths μ0. x-axis, Difference between μ0 and μ*0, which is the bifurcation point at which the saddle steady state appears and whose value depends on the AMPA:NMDA ratio. C and D have the same symbolic notations as in B. Error bars indicate SD.
For an AMPA:NMDA ratio larger than 50:50, we found that the network became dynamically unstable leading to large-amplitude oscillations (supplemental Fig. 4, supplementary information F, available at www.jneurosci.org as supplemental material). This is because, when recurrent excitation is dominated by the AMPA receptors, the interplay between fast positive feedback (τAMPA = 2 ms) and slower negative feedback (mediated by the GABAA synapses, τGABAA = 5 ms) naturally gives rise to oscillatory instability (Wang, 1999; Compte et al., 2000; Tegner et al., 2002; Brunel and Wang, 2003). To assess time integration solely mediated by AMPA receptors, we considered the artificial case where τGABA = τAMPA = 2 ms. (To preserve the same steady-state behavior, the term νI in Equation 7 is multiplied by a factor of τGABA/τAMPA. To compensate for the loss in NMDA synaptic currents, the peak synaptic conductances for AMPA, grec, AMPA, between excitatory cells and excitatory-to-inhibitory connections were changed to 0.000237 and 0.000289 μS, respectively. The recurrent strength w+ was also raised to 2.4.)
Figure 8 displays simulations with 100% AMPA at recurrent synapses. We found that because of the very fast dynamics (Fig. 8A), it is virtually impossible to achieve an integration time of >200 ms, even when the parameter μ0 was fine-tuned to be very close to the bifurcation point μ*0 (Fig. 8B).
Decision time with only AMPA at recurrent synapses (c′ = 0%). A, A sample time course with very fast dynamics (t = 0 is stimulus onset). μ0 = μ*0 + 15 Hz, where μ*0 is the value of μ0 when the saddle point is just created. Bottom horizontal bar denotes the duration where the firing rates are indistinguishable; trajectory lies near the stable manifold of the saddle point. B, Reaction time as a function of μ0 − μ*0. Note that reaction times longer than 200 ms can hardly be realized even with parameter fine-tuning. Error bars indicate SD.
These results show that the NMDA receptors at recurrent synapses are important to slow time integration in the model.
Decision making and working memory
Time integration with a mnemonic delay period
We have mentioned previously that, like the original model (Wang, 2002), our reduced model can also simulate the delayed response version of the visual motion discrimination task (Shadlen and Newsome, 2001) in which the monkey is required to withhold the response across a delay period of a few seconds when the choice must be maintained actively in working memory. As shown in Figure 9, the model neuronal activities are comparable with the observations from LIP cells (Shadlen and Newsome, 2001). In particular, in response to a 1 s motion stimulus presentation, the ramping of the neural activity is faster at higher coherence levels (Fig. 9A). After removal of stimulus, a sustained persistent activity of the decision is stored in working memory. Neurons whose response fields are opposite to that of the saccadic target eventually show a ramping down activity. This is expected because of the effective mutual inhibition between competing neural groups. Note that, similar to the reaction time simulations, the time courses of the two neural populations are initially indistinguishable for a hundred of milliseconds, before they eventually split apart.
Neural firing activity in delayed-response task in which a motion stimulus presentation is followed by a memory period. A, Sample time courses within a trial with different coherence levels. Black and dashed lines are for saccades moving toward and away from the response field of the neuron, respectively. The coherence level is shown at the top of each panel. Shaded regions, Motion viewing period. The black horizontal line at the bottom indicates the time epoch when the two firing rates are indistinguishable. B, Dependence of neural activity on motion strength in different epochs. Opened and filled circles, Saccades toward and away from the response field of the neuron. The largest dependence on the motion strength, as well as the greatest difference in the two neural responses, correspond to the late phase (0.5–1 s epoch) of stimulus presentation. In the third epoch (early delay period), there is a still a residual effect of the dependence of neural response on motion strength. Figures are calculated using correct trials only and averaged over 2000 trials.
In Figure 9B, we show the neural response for each coherence level averaged over a duration of 0.5 s. Similar to the physiological observations [Shadlen and Newsome (2001), their Fig. 9], the difference of the firing rates between the two choices is small in the initial epoch, and becomes the largest later, during stimulus presentation. It is also the second time epoch, 0.5–1 s after motion onset, that shows the greatest sensitivity of the neural response to the coherence level. During the delay period, the input is withdrawn and the system is approaching a steady state of persistent activity; hence there is a decrease in the difference in neural responses with respect to coherence levels. Interestingly, like LIP neurons, the firing activity of our model during the early epoch (1–1.5 s) of the delay period still shows a slight residual dependence on motion strength (Fig. 9B), which disappears later (1.5–2 s) when the mnemonic steady state is reached.
Effect of external input on network behavior
We have seen that the model system has different steady states of high neural activities in the presence or absence of an external stimulus. To gain intuition about how these two steady states are related, we computed the steady states of the system as a function of a stimulus parameter, such as the amplitude of an unbiased input μ0 (Fig. 10). In this “bifurcation diagram, ” the middle curve corresponds to symmetric steady states (when S1 = S2), whereas the upper and lower branches correspond to asymmetric states (when S1 is high and S2 is low, or vice versa). Solid and dashed lines denote stable and unstable (saddle-like) steady states, respectively. It is immediately clear that the steady states do not simply vary quantitatively and gradually. Instead, they can change their stability, disappear, or emerge, as μ0 is varied. Briefly, the network goes from being monostable to bistable (coexistence of a resting state and mnemonic persistent states) as μ0 is increased (Fig. 10). Additional increase of μ0 to a critical value of μ0* ≈ 7 Hz causes the spontaneous steady state to lose stability and becomes a saddle point, a bifurcation of the system.
Bifurcation diagram of a selective population with stimulus strength μ0 as parameter (c′ = 0%). Bold lines, Stable steady states; dashed lines, saddle steady states. Spontaneous state before stimulus presentation is denoted by the filled square. With a μ0 = 30 Hz stimulus, the spontaneous stable state loses stability, and a saddle steady state appears (open square). The state either goes toward the upper or lower stable state (filled circles). The population wins the competition if the upper branch is chosen, and loses otherwise. When stimulus is removed, hysteresis of the upper stable branch allows the activity to persist (memory storage of a decision choice). Arrow with an asterisk, Point where spontaneous state loses stability. Arrow with double asterisks, Saddle point turns into an attractor.
At μ0** ≈ 47 Hz, another bifurcation occurs where the saddle point changes stability and is changed to a symmetric stable steady state with high S. In contrast, the asymmetric attractor states exist only for a finite range of μ0 values: a sufficiently large unbiased input μ0, either negative or positive, would lead to the disappearance of the asymmetric attractor dynamics, hence the loss of the ability of the system to compute a categorical choice through winner-take-all competition.
The model dynamics that was shown in Figure 4, A and B, can be alternatively viewed as follows. At rest, the system is in the symmetric state (filled square) at μ0 = 0. When a sufficiently large external stimulus is presented, the system is suddenly brought to another point in the bifurcation diagram (e.g., at μ0 = 30 Hz), where this symmetric steady state is unstable (open square). Because of this instability, the population activity is attracted to either the upper or lower stable branch. If the system goes to the upper branch, then this population “wins” while the other population “loses” (filled circles). Furthermore, when the external stimulus is removed from the network, the state of either population stays at the same upper or lower stable branch (filled triangles) because the two branches still exist at μ0 = 0. This hysteresis phenomenon acts as a short-term memory of the decision/choice.
Two distinct regimes of decision-making dynamics
We have shown that, with a strong self-excitation within the network, both decision making and working memory can be achieved. In particular, we saw that they are related to the same “branch” of stable steady states in a bifurcation diagram. We will now gradually decrease this self-excitatory component, in the form of w+, and study how time integration is affected by it.
To some degree, we expect the effect of varying the recurrent strength w+ would be similar to that with varying the stimulus μ0. Given the same stimulus μ0, the lower w+ is, the slower the ramping neural activity, and the longer the decision time (Fig. 11A). This can be explained by the fact that both w+ and μ0 have a direct influence on the two competing neural populations. Thus, lowering either w+ or μ0 causes the saddle point to approach the spontaneous state, slowing the trajectories and thus the integration time (Fig. 11B). We also assessed quantitatively how w+ affects the speed and accuracy of decisions over a range of c′. We found that, with a larger w+, the decision became faster across all values of c′ (Fig. 11C). Note that for large c′, the decision time is less sensitive to w+, as the system is strongly attracted toward the correct attractor. At the same time, the performance was worse (Fig. 11D) with a larger w+. Quantitatively, when w+ was increased (1.68, 1.7, or 1.72), the threshold of the psychometric function (α) increased (4.10, 5.50, or 8.52% coherence) while the slope (β) was 1.31, 1.25, and 1.36.
Dependence of integration time on the relative strength of recurrent excitation. A, Sample time courses: faster ramping activity with stronger recurrent strengths, w+. B, The unstable time constant of a saddle point dominates the dynamics when recurrent strength w+ is weak. C, Reaction time decreases with increasing w+. Error bars indicate SD. D, Accuracy of performance decreases with increasing w+. Data are fit with a Weibull function.
We have demonstrated how the reaction time and performance can be affected by the recurrent strength w+ and, earlier on, the stimulus input strength μ0. To further study how the combined interaction between w+ and μ0 can affect time integration, we constructed a state diagram of the reduced decision-making model in the space of both parameters w+ and μ0 (Fig. 12). We found that this parameter plane is divided into three regions: monostability (with only one symmetric steady state), competition (with only asymmetric attractor states), and bistability (coexistence between one symmetric and two asymmetric attractor states). Depending on the recurrent wiring property (w+ value), the system behaves as function of a stimulus μ0 in four different ways. On one hand, in a network with weak recurrent connections (w+ < 1.58), the network cannot serve working memory because of the lack of bistability, nor can it perform decision-making computation, because no input can bring it into the competition region of the state diagram (Fig. 12, inset, regime I). On the other hand, in a network with sufficiently strong recurrent connections (1.6 < w+ < 1.86), the network behaves the same way as with the standard set of parameters that we have observed before; by increasing μ0 (to within a suitable range), we can bring the network from a spontaneous state in the bistable (left red) region to a competitive (blue) region where only two competing choice attractors exist. And when stimulus is removed, hysteresis in the same bistable region, allows the decision to be stored in short-term memory (Fig. 12, inset, regime III, or Fig. 10). Hence, the same network can subserve both decision-making computation and working memory. In the extreme case when w+ is very high (w+ > 1.86) (Fig. 12, inset, regime IV), there is no competition in the sense that, for any μ0, there is always a stable symmetric state that coexists with the asymmetric attractors. This regime will not be considered further in the present study.
Distinct modes of operation in the two parameter space with zero coherence. In general, there are three types of regions. Bistable (red) region, A symmetric and two asymmetric attractors coexist; blue competition region, one saddle with two asymmetric competing attractors; monostable region, only one attractor. Depending on the strength of recurrent excitation w+, the network responds to a stimulus (of suitable intensity μ0) in four different ways, shown as regimes I, II, III, IV in insets. Regime I and II do not support working memory (of decision). Regime I, No decision making nor memory. Regime II, The network can produce a binary decision during stimulation but cannot store it in working memory. Regime III, The network is capable of both decision-making computation and working memory (our standard parameter set). Regime IV, For any μ0, there is always a stable symmetric stable state. Dark and dashed branches denote loci of stable and unstable steady states, respectively. A and AS are labels for branches with symmetric and asymmetric steady states, respectively.
In a network with moderate recurrent connections (1.58 < w+ < 1.6), with a suitable range of μ0, the network can also be brought to a region of competition conducive for decision making. An example is shown in Figure 13. However, on the removal of stimulus, the network lies outside the bistable region and thus does not display hysteresis (Fig. 12, inset, regime II). As a result, the network is unable to sustain any persistent activity without stimulus (Fig. 13A, B), and there is no memory storage of the decision.
Decision making without short-term memory in a network with low recurrent strength (w+ = 1.59). A, A typical trial showing slow ramping up activity for both neural populations. Note that the two firing rates are indistinguishable for many hundreds of milliseconds (indicated by the black horizontal bar), before they separate by a modest difference. Stimulus is applied from time 0–2 s (gray region). B, Phase-plane without stimulus has only one low stable attractor. C, The response of the system to a stimulus with zero coherence, in two trials, plotted as trajectories in the decision space (red and blue). Stimulus intensity is μ0 = 45 Hz lasting for 2 s. D, Comparison between stable and unstable time constants of the saddle-type unstable steady state (μ0 = 45 Hz).
Thus, our model in regime II can perform decision computation (Fig. 13A, C) without working memory function, similar to previously studied connectionists models (Brown and Holmes, 2001; Usher and McClelland, 2001). Furthermore, we found that, for the whole range of μ0, the unstable time constant of the saddle point is of a few seconds, an order of magnitude greater than the stable time constant (Fig. 13D). Time integration is then dominated by the unstable manifold (Fig. 13B). Thus, qualitatively, the network dynamics may be approximated by a picture in which the system initially collapses onto the unstable manifold, followed by a one-dimensional diffusion-like dynamics for the difference S1 − S2 (Brown and Holmes, 2001). However, quantitatively, the dynamics along the stable manifold is still slow (a few hundreds of milliseconds) (Fig. 13A, D), and thus cannot be neglected. Moreover, with the parameter values explored, the activities are very low (Figs. 12, inset II; 13A), and reaction times are too long (Fig. 13A) in comparison with the monkey behavioral data. Finally, it is worth noting that regime II is realized in our model with a relatively small range of w+ values (between 1.58 and 1.6). For these reasons, we conclude that, although our model in regime II can mathematically be reduced to a diffusion-like model, it does not quantitatively reproduce the reaction times and firing rates observed in the experiment of Roitman and Shadlen (2002).
Discussion
A dynamical system approach to decision making
To elucidate the wiring properties and neural dynamics of a cortical microcircuit subserving decision making, we developed a two-variable model that was derived from a biophysically based cortical network model of decision making, consisting of thousands of spiking neurons. This was done using a mean-field approach and through a number of approximations. Such a reduction represents a significant step in bridging the gap between the biologically based model of Wang (2002) and abstract mathematical models, and allow for a comparison between our model and previously proposed models of decision making. Like other models (Usher and Cohen, 1999; Brown and Holmes, 2001; Usher and McClelland, 2001; Machens et al., 2005), our model is a nonlinear dynamical system that can be conceptualized as a circuit of neural populations coupled effectively by mutual inhibition. However, it is important to emphasize that recurrent excitation plays an indispensible role in producing reverberatory dynamics underlying winner-take-all competition in our model. Furthermore, the recurrent excitation is dominated by a slow (and saturating) recurrent synaptic dynamics (of NMDA receptors). Reciprocal inhibition between choice-selective neural pools is mediated by feedback inhibition in the local circuit, rather than through cross-inhibition along feedforward input pathway as assumed in a previous model (Shadlen and Newsome, 2001; Mazurek et al., 2003). A prediction from the local inhibition architecture, but not by the feedforward model, is that, if an excitatory perturbation is applied to neural pool 1 during sensory stimulation, not only will it accelerate the decision time for choice 1, but also slow down the decision time for choice 2. This was indeed found to be the case in a recent study where microstimulation was applied to LIP neurons in the reaction time task of visual motion discrimination (Hanks and Shadlen, 2004).
Our reduced model not only reproduces salient experimental findings of Roitman and Shadlen (2002) and simulations of the original model of Wang (2002), it also provides a neurodynamical framework for a deeper understanding of time integration and categorical choice computation in the decision process. In particular, the phase-plane analysis sheds insights into how the system works in several ways. First, it shows how, in response to a stimulus, the network changes its configuration into a mode for winner-take-all competition. Second, it revealed that local dynamics in the vicinity of a saddle point plays a key role in controlling the slow temporal course of stimulus integration, especially when the stimulus is weak and the coherence is low or zero. Third, it offers a precise explanation as to why decision times are slower in error trials than in correct trials, as in the monkey experiment (Roitman and Shadlen, 2002). Fourth, we found that when the coherence is above a critical value (∼70%), one of the competing attractor states disappears, so that the system is forced to go to the other attractor. Interestingly, in simulations of the original model in which the sign of the input was reversed during stimulation (Wang, 2002), it was found that a decision made by a stimulus presentation can always be reversed by a second stimulus with opposite sign, if the reversing signal strength was above c′ = 70%, no matter how late the reversing signal was applied after stimulus onset. Here, this observation is nicely explained by the fact that, in the presence of a stimulus with sufficiently strong motion strength c′, only one attractor exists in the phase space so that the choice by the system is unique in all trials.
Slow recurrent dynamics underpinning integration time
An important result of this work is that when recurrent excitation has a significant NMDA component, temporal integration is generally slow, leading to accurate decisions. We asked ourselves whether this is necessarily the case, by gradually substituting NMDARs by faster AMPARs at recurrent synapses. Consistent with previous observations in such a network model (Wang, 1999; Compte et al., 2000; Tegner et al., 2002), we found that the network becomes highly unstable when the local reverberation is largely mediated by AMPARs. Moreover, when we artificially prevented oscillatory instability (under the assumption that AMPA current and GABAA receptor-mediated inhibitory current have a similar decay time constant), the model cannot reproduce long decision times comparable with behavioral data, even with fine-tuning of parameters to maximize the integration time of the system. This results suggests that NMDARs at recurrent synapses in a decision circuit may be critically important. Experimentally, it would be interesting to see whether motion discrimination becomes more impulsive and less accurate when NMDAR antagonists are applied to LIP in behaving nonhuman primates.
Future work is needed to examine quantitative properties of the NMDAR-mediated synaptic current, such as the temperature dependence of the current kinetics (Hestrin et al., 1990) or subunit composition of NMDARs (McBain and Mayer, 1994), and their influence on time integration of a cortical circuit. Moreover, other slow cellular processes, like calcium-dependent inward current (Egorov et al., 2002; Durstewitz, 2003; Major and Tank, 2004) or short-term synaptic plasticity (Abbott and Regehr, 2004) may contribute to time integration as well as persistent activity. It would be interesting to incorporate these biophysical processes in a decision-making network model.
Neuronal versus diffusion-like models
Previous work has shown that connectionist models with mutually inhibitory neural pools can be reduced to a one-dimensional diffusion model, under two conditions. First, the system can be approximately described by a single dynamical variable, the firing rate difference between the two neural populations. Second, the time constant for this dynamical variable is arbitrarily long, which typically requires fine-tuning of parameters (Brown and Holmes, 2001; Usher and McClelland, 2001). Our biophysically based spiking neuron network model provides an opportunity to assess whether, or how, these conditions can be satisfied in a realistic cortical circuit. We found that the first condition breaks down into two requirements: (1) that time integration is dominated by dynamics near an unstable saddle point and (2) that the stable time constant of the saddle is negligible.
For our model, we found that (1) holds true for stimuli with very low or zero coherence c′, which are the most interesting situations where time integration is crucial. In general, with a reasonably strong stimulus intensity (μ0), the stable time constant (τstable) is larger than the unstable one (τunstable) of the saddle. This is reflected by a biphasic time course of neural Apopulation firing rates in response to stimulation. With weak μ0, τunstable becomes longer than τstable, but even in that case τstable is ∼200–300 ms, so quantitatively its contribution to integration time cannot be ignored. Therefore, (2) is generally not satisfied in our model. As to the second condition, fine-tuning amounts to adjusting the system to a bifurcation point, where a time constant (τunstable in our model) diverges to infinity. However, we found that this is unnecessary. Our model endowed with slow reverberation (primarily mediated by the NMDARs) can reproduce the decision times of many hundreds of milliseconds up to a second robustly, without fine-tuning of parameters.
One may argue that the condition (2) for τstable to be negligibly short could be fulfilled when the recurrent excitation is dominated by AMPARs rather than NMDARs (see supplementary information E, available at www.jneurosci.org as supplemental material) (Figs. 7, 8A). If so, then τunstable could indeed dominate the network dynamics. However, because in this scenario the system does not have intrinsically slow time constants, long integration time is impossible without fine-tuning of parameters. In fact, despite our efforts, we were hardly able to realize reaction time of >200 ms even when the system is extremely fine-tuned to a bifurcation point. For these reasons, we are in favor of a two-dimensional dynamical system model that robustly performs NMDAR-dependent time integration.
Our model naturally explains the observation that longer reaction times in error trials compared with correct trials in the monkey experiment of Roitman and Shadlen (2002) as well as in most human reaction time studies (Ratcliff, 1978; Luce, 1986), a feature that cannot be reproduced by the diffusion model unless some additional ingredients are added to it (Ratcliff and Rouder, 1998; Mazurek et al., 2003). Neurophysiologically, there is evidence that LIP neurons recorded from behaving monkeys also exhibit biphasic temporal dynamics similar to our model (Shadlen and Newsome, 2001; Huk and Shadlen, 2005). For a comparison, the two neural population activities of our model would correspond in physiological data to spike activities of LIP cells pooled over trials according to the animal’s choice (toward or away from the response field of the cell). Qualitatively consistent with our model, after the stimulus onset these two time courses are initially indistinguishable (for 100 ms or more) before diverging from each other (one increases while the other declines over time). This observation remains preliminary, and it would be desirable to analyze it in more detail. It is hoped that additional experiments will be performed to further test our two-dimensional neuronal model versus the one-dimensional model both behaviorally and neurophysiologically.
Distinct modus operandi of the model
We mapped out a two-parameter (μ0, w+) state diagram of the decision-making model and found that there is a finite parameter region suitable for competition, where a categorical decision must occur because the only stable attractors are the choice attractors. This range of w+ values is finite and preliminary observations (data not shown) indicate that it can be enlarged with a higher AMPA:NMDA ratio at the recurrent synapses.
As long as w+ is not too small, a suitable stimulus (μ0 ≠ 0) can bring the network into the competition region. This is the case whether our network at rest (with μ0 = 0) is either monostable (regime II) or bistable (regime III), depending on the strength of w+. In regime II, our network is capable of decision computation without working memory, but the corresponding range of w+ values is small. Also, in this regime neuronal firing rates are fairly low compared with LIP data, with the parameters we have used. A possible way of increasing the firing rates is by having a higher proportion of AMPA at the recurrent connections. It remains to be determined, in future work, whether this regime can be realized in a more robust manner, for example with parameter changes or by adding new ingredients into the model. Experimentally, it is an important and still open question of whether the LIP local circuit is capable of sustained persistent activity by itself, and operates as an attractor network.
Wang (2002) and the present study show that a single local network can naturally serve both working memory and decision making. The same conclusion was reached by Machens et al. (2005), in a modeling work concerned with a delayed somatosensory discrimination experiment (Romo et al., 1999). In that task, two stimuli presented separately in time must be compared, for the animal to reach a categorical perception. Thus, there is no explicit need for slow time integration of sensory stimulus. In contrast, the first stimulus must be held in short-term memory across the delay, as an analog quantity, which is encoded monotonically and stored by a line attractor in the model (Miller et al., 2003; Machens et al., 2005). This is different from the visual discrimination experiment of Shadlen and Newsome (2001), and our model, in which a sensory information is accumulated over time and a categorical decision is reached before a mnemonic delay period. Therefore, whereas the integrated information is an analog quantity, the stored information is discrete (binary choice). It is conceivable that a different design of the visual motion discrimination task may require integration of analog sensory data across “temporal delays, ” which could be substantiated by a continuous attractor network.
Our results, in consonance with those of Machens et al. (2005), show that an attractor dynamical system is readily re-configurable by external inputs. So far, we have viewed μ0 as representing sensory information conveyed from the MT neurons to the LIP neurons. The value of μ0 may depend on contrast or speed of motion stimulus. But because μ0 is a generic input to the LIP network, we need not restrict the interpretation of it to bottom–up input. We can also think of μ0 as partly originating from top–down modulatory pathways. For instance, “executive control” signals from the prefrontal cortex can bring an LIP network in and out of the bistability regime, so that the system can operate as an attractor network, or not, depending on behavioral demands.
Future directions
The attractor model in Wang (2002) and the present reduced model are incomplete in several aspects. First of all, instead of having ∼2–3 Hz spontaneous firing rate in the models, the firing rates before stimulus presentation in Shadlen and Newsome (2001) and Roitman and Shadlen (2002) are much higher (>15 Hz). This might be attributable to the presence of the two target stimuli that signal the two alternative choices to the monkey in the experiment. We have tested, and confirmed, this idea by adding a strong input that is symmetrical to both neural populations, before a motion stimulus is presented. Detailed results are beyond the scope of this paper and will be reported elsewhere.
Second, both our model and the diffusion-type model assumes that the event of threshold crossing in a decision network, can be read out by a downstream system leading to a motor response. How this is instantiated in the brain remains unknown [for a possible biological mechanism, see Lo and Wang (2004)].
Third, our model does not capture effects of previous trials on the performance of the present trial (Brown and Holmes, 2001; Bogacz et al., 2003), which may be caused by priming (Carpenter and Williams, 1995; Fecteau and Munoz, 1999; Yang and Shadlen, 2004, 2005) or reflect an optimization of overall rewards (Glimcher, 2003; Brown et al., 2005; Sugrue et al., 2005). The cellular mechanisms underlying trial-by-trial choice correlations presumably involve learning at the synaptic level (Fusi et al., 2005), which should thus be incorporated and explored in a more complete neural network model of perceptual decision making in reaction time tasks.
Appendix
Reduced two-variable model without AMPA at recurrent synapses
The reduced two-variable model is in its simplest form if we assume that there is no AMPA at the recurrent synapses. In this case, the system can be completely described by the following equations: 
where i = 1, 2 labels the selective population. Parameter values for the input–output function are a = 270(VnC)−1, b = 108 Hz, and d = 0.154 s. The kinetic parameters are γ = 0.641, τS = 100 ms, and τAMPA = 2 ms. The synaptic couplings are JN, 11 = JN, 22 = 0.2609 nA, JN, 12 = JN, 21 = 0.0497 nA, and JA, ext = 5.2 × 10−4 nA · Hz−1. The overall effective external input is I0 = 0.3255 nA, noise amplitude is σnoise = 0.02 nA, and the stimulus is μ0 = 30 Hz. The value of σnoise is chosen such that the psychometric and chronometric functions are close to that of Roitman and Shadlen (2002). A Matlab code for simulation and an XPPAUT code for phase-plane analysis can be obtained from the authors on request.
Footnotes
-
This work was supported by National Institutes of Health Grants MH062349 and DA016455 and by the Swartz Foundation. We are grateful to Stefano Fusi and Alireza Soltani for helpful discussions. We also thank D. J. Barraclough, N. Brunel, E. S. Carter, F. Liu, and P. Miller for comments on this manuscript.
- Correspondence should be addressed to Xiao-Jing Wang at the above address. Email: xjwang{at}brandeis.edu
