Abstract
Recent physiological studies using behaving monkeys revealed that, in a twoalternative forcedchoice 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 twovariable 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 winnertakeall 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 twovariable neural model offers a simple yet biophysically plausible framework for studying perceptual decision making in general.
 sensory discrimination
 intraparietal cortex
 reaction time
 computational modeling
 dynamical systems
 NMDA
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 ramptothreshold 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 onedimensional 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 finetuning of parameters cancels out the leakage by network recurrent dynamics (Brown and Holmes, 2001; Usher and McClelland, 2001; Bogacz et al., 2003).
Although the diffusiontype 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 winnertakeall 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 meanfield 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 meanfield 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 “meanfield” approach to reduce the model. This approach has been used to analytically study spiking neuronal network models comprised of integrateandfire 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 selfconsistently. 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 alltoall 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.
Moreover, the firing rate r of a leaky integrateandfire (LIF) neuron receiving a noisy input current is given by the firstpassage time formula (Ricciardi, 1977; Amit and Tsodyks, 1991; Renart et al., 2003): (1) where φ is a function of the total synaptic input current I_{syn}. τ_{m} is the membrane time constant. V_{th} is the spiking threshold for the membrane voltage, V_{reset} is the reset voltage, τ_{ref} is the refractory period, σ_{V} is the membrane potential SD, and V_{ss} = V_{L} + I_{syn}/g_{L}. 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. I_{syn} is the total synaptic input to a single cell, and c_{E, I} is the gain factor. g_{E, I} is a noise factor that determines the shape of the “curvature” of φ. If g_{E, I} is large, φ would act like a linearthreshold function with I_{E, 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 firstpassage time formula Equation 1 of a singlecell LIF model driven by AMPA receptormediated 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, I_{E} = 125 Hz, g_{E} = 0.16 s, and c_{E} = 310(VnC)^{−1}; and for interneurons, I_{I} = 177 Hz, g_{I} = 0.087 s, and c_{I} = 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 meanfield 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. r_{i}(t) is the instantaneous mean firing rate of the presynaptic excitatory population i, and r_{I}(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 S_{i}.
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 deltalike 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].
Phaseplane reduction
The model can be further reduced to a twovariable 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 singlecell 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 g_{2} = 2 and r_{0} = 11.5 Hz.
With the linearity in the response of the inhibitory population, the selfinhibitory term can now be easily absorbed into the factors g_{2} and r_{0}. In particular, if we define J_{II} as the selfinhibitory synaptic coupling, then selfinhibition is expressed as a term −J_{II}φ in I_{syn}. Grouping the two φ dependent terms together, it is clear that selfinhibition term effectively lowers the mean firing rate of interneurons, r_{I}, by a factor of 1 + (c_{I}/g_{2})J_{II}. This helps to simplify the calculations by removing selfconsistency 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 GABA_{A} 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 S_{NMDA}, 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 GABA_{A} gating variables reach their steady states much faster than that of NMDA receptors, which means that the average gating variables of AMPA and GABA_{A} 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 twovariable 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 S_{NMDA} and τ_{S} for τ_{NMDA}. The firing rates r_{1} and r_{2} are given by Equation 2: (12) (13) (14) (15) where I_{i} represents the visual motion stimulus to the population i and depends on the motion strength (see Results). I_{noise, i} is a noise term, and I_{0} is the mean effective external input common to both populations. Because of the background input into nonselective excitatory cells and interneurons, I_{0} includes not only direct background input to a selective population but also indirect background inputs from these nonselective cells. I_{syn, i} is the total synaptic current from both recurrent connections and inputs fed from outside the local network. The coefficients J_{N, ij} and J_{A, ij} are effective coupling constants from neuron j to i mediated by NMDAR and AMPAR, respectively. The negative sign in front of J_{N, 12}, J_{N, 21}, J_{A, 12}, and J_{A, 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 J_{N, I→E}J_{E→I}(S_{1} + S_{2}) and J_{A, I→E}J_{E→I}(r_{1} + r_{2}). Thus, for example, in I_{syn, 1}, the S_{2}dependent term is of the form (J_{N, E→E} − J_{N, I→E}J_{E→I})S_{2} ≡ −J_{N, 12}S_{2} with J_{N, 12} = J_{N, I→E}J_{E→I} − J_{N, E→E}· J_{N, 21}, J_{A, 12}, and J_{A, 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: r_{i} = φ(I_{syn, i}), where I_{syn, i} depends on both r_{1} and r_{2}. This problem was solved by finding an effective singlecell 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 selfcontained twovariable system for S_{1} and S_{2}: (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 twovariable model is as follows: J_{N, 11} = 0.1561 nA = J_{N, 22}, J_{N, 12} = 0.0264 nA = J_{N, 21}, J_{A, 11} = 9.9026 × 10^{−4} nC = J_{A, 22}, J_{A, 12} = 6.5177 × 10^{−5} nA · Hz^{−1} = J_{A, 21} and I_{0} = 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 meanfield approach does not include timevarying noise that plays a critical role in the spiking neural network. To amend this, we added a noise term I_{noise} 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 σ_{noise}^{2} 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.
Phaseplane 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 trialbytrial 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 J_{A, 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 twovariable model replicates fairly well the psychometric function (behavioral performance) and chronometric function (reaction time of correct and error trials) of both the original largescale 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.
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 phaseplane analysis of the model system. This is done by first setting the dynamical equations dS_{1}/dt = G_{1}(S_{1}, S_{2}) = 0 and dS_{2}/dt = G_{2}(S_{1}, S_{2}) = 0, and then plotting these two lines in the (S_{1}, S_{2}) 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 S_{i} of population i is a monotonic, increasing function of the population firing rate r_{i} and their steady states are related by the simple monotonic Equation 8 (see Materials and Methods), we expect the nullclines in the (r_{1}, r_{2}) 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 lowerleft 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 selfsustained 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 stimulusselective persistent activity states, which can be internally maintained across a delay period.
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 saddletype 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 twoalternative forcedchoice 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 phasespace 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 decisionmaking dynamics in terms of a onedimensional “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, S_{1} − S_{2} or r_{1} − r_{2}). 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 shortterm 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 onedimensional, although we have seen that the dynamics of the model should be understood in the twodimensional phase space. Below, we will address the question of whether the system can be reduced to a onedimensional 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.
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 − S_{saddle}, between S = (S_{1}, S_{2}) and the saddle S_{saddle} can be written as follows: (22) where v_{1} and v_{2} 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 a_{1} and a_{2} 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 a_{2} = 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, a_{1} = 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(t_{stay}/τ_{unstable}), which yields t_{stay} ≃ τ_{unstable}log(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 onedimensional 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 S_{1} − S_{2} (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.
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 onedimensional 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 finetuning 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 twovariable model, which does not take into account the fast AMPA dynamics. Instead, we investigated the full, 11variable, 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, g_{AMPA} and g_{NMDA}, 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).
For an AMPA:NMDA ratio larger than 50:50, we found that the network became dynamically unstable leading to largeamplitude 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 GABA_{A} 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 steadystate 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, g_{rec, AMPA}, between excitatory cells and excitatorytoinhibitory 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 finetuned to be very close to the bifurcation point μ*_{0} (Fig. 8B).
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.
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 S_{1} = S_{2}), whereas the upper and lower branches correspond to asymmetric states (when S_{1} is high and S_{2} is low, or vice versa). Solid and dashed lines denote stable and unstable (saddlelike) 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.
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 winnertakeall 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 shortterm memory of the decision/choice.
Two distinct regimes of decisionmaking dynamics
We have shown that, with a strong selfexcitation 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 selfexcitatory 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.
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 decisionmaking 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 decisionmaking 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 shortterm memory (Fig. 12, inset, regime III, or Fig. 10). Hence, the same network can subserve both decisionmaking 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.
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.
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 onedimensional diffusionlike dynamics for the difference S_{1} − S_{2} (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 diffusionlike 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 twovariable 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 meanfield 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 winnertakeall competition in our model. Furthermore, the recurrent excitation is dominated by a slow (and saturating) recurrent synaptic dynamics (of NMDA receptors). Reciprocal inhibition between choiceselective neural pools is mediated by feedback inhibition in the local circuit, rather than through crossinhibition 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 phaseplane 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 winnertakeall 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 GABA_{A} receptormediated inhibitory current have a similar decay time constant), the model cannot reproduce long decision times comparable with behavioral data, even with finetuning 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 NMDARmediated 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 calciumdependent inward current (Egorov et al., 2002; Durstewitz, 2003; Major and Tank, 2004) or shortterm 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 decisionmaking network model.
Neuronal versus diffusionlike models
Previous work has shown that connectionist models with mutually inhibitory neural pools can be reduced to a onedimensional 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 finetuning 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, finetuning 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 finetuning 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 finetuning of parameters. In fact, despite our efforts, we were hardly able to realize reaction time of >200 ms even when the system is extremely finetuned to a bifurcation point. For these reasons, we are in favor of a twodimensional dynamical system model that robustly performs NMDARdependent 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 twodimensional neuronal model versus the onedimensional model both behaviorally and neurophysiologically.
Distinct modus operandi of the model
We mapped out a twoparameter (μ_{0}, w_{+}) state diagram of the decisionmaking 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 shortterm 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 reconfigurable 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 diffusiontype 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 trialbytrial 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 twovariable model without AMPA at recurrent synapses
The reduced twovariable 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 J_{N, 11} = J_{N, 22} = 0.2609 nA, J_{N, 12} = J_{N, 21} = 0.0497 nA, and J_{A, ext} = 5.2 × 10^{−4} nA · Hz^{−1}. The overall effective external input is I_{0} = 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 phaseplane analysis can be obtained from the authors on request.
Footnotes
 Received September 6, 2005.
 Revision received December 8, 2005.
 Accepted December 11, 2005.

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 XiaoJing Wang at the above address. Email: xjwang{at}brandeis.edu