Drift–diffusion models for multiple-alternative forced-choice decision making

The canonical computational model for the cognitive process underlying two-alternative forced-choice decision making is the so-called drift–diffusion model (DDM). In this model, a decision variable keeps track of the integrated difference in sensory evidence for two competing alternatives. Here I extend the notion of a drift–diffusion process to multiple alternatives. The competition between n alternatives takes place in a linear subspace of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$n-1$\end{document}n−1 dimensions; that is, there are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$n-1$\end{document}n−1 decision variables, which are coupled through correlated noise sources. I derive the multiple-alternative DDM starting from a system of coupled, linear firing rate equations. I also show that a Bayesian sequential probability ratio test for multiple alternatives is, in fact, equivalent to these same linear DDMs, but with time-varying thresholds. If the original neuronal system is nonlinear, one can once again derive a model describing a lower-dimensional diffusion process. The dynamics of the nonlinear DDM can be recast as the motion of a particle on a potential, the general form of which is given analytically for an arbitrary number of alternatives.


Introduction
Perceptual decision-making tasks require a subject to make a categorical decision based on noisy or ambiguous sensory evidence. A computationally advantageous strategy in doing so is to integrate the sensory evidence in time, thereby improving the signal-to-noise ratio. Indeed, when faced with two possible alternatives, accumulating the difference in evidence for the two alternatives until a fixed threshold is reached is an optimal strategy, in that it minimizes the mean reaction time for a desired level of performance. This is the computation carried out by the sequential probability ratio test devised by Wald [1], and its continuous-time variant, the drift-diffusion model (DDM) [2]. It would be hard to overstate the success of these models in fitting psychophysical data from both animals and human subjects in a wide array of tasks, e.g. [2][3][4][5], suggesting that brain circuits can implement a computation analogous to the DDM.
At the same time, neuroscientists have characterized the neuronal activity in cortical areas of monkey, which appear to reflect an integration process during DM tasks [6], although see [7]. The relevant computational building blocks, as revealed from decades of in-vivo electrophysiology, seem to be neurons, the activity of which selectively increases with increasing likelihood for a given upcoming choice. Attractor network models, built on this principle of competing, selective neuronal populations, generate realistic performance and reaction times; they also provide a neuronal description which captures some salient qualitative features of the in-vivo data [8,9].
While focus in the neuroscience community has been almost exclusively on twoalternative DM (although see [10]), from a computational perspective there does not seem to be any qualitative difference between two or more alternatives. In fact, in a model, increasing the number of alternatives is as trivial as adding another neuronal population to the competition. On the other hand, how to add an alternative to the DDM framework does not seem, on the face of things, obvious. Several groups have sought to link the DDM and attractor networks for two-alternative DM. When the attractor network is assumed linear, one can easily derive an equation for a decision variable, representing the difference in the activities of the two competing populations, which precisely obeys a DDM [11]. For three-alternative DM previous work has shown that a 2D diffusion process can be defined by taking appropriate linear combinations of the three input streams [12]. The general nalternative case for leaky accumulators has also been treated [13]. In the first section of the paper I will summarize and build upon this previous work to illustrate how one can obtain an equivalent DDM, starting from a set of linear firing rate equations which compete through global inhibitory feedback. The relevant decision variables are combinations of the activity of the neuronal populations, and which represent distinct modes of competition. Specifically, I will propose a set of "competition" basis functions which allow for a simple, systematic derivation of the DDMs for any n. I will also show how a Bayesian implementation of the the multiple sequential probability ratio test (MSPRT) [14][15][16] is equivalent in the continuum limit to these same DDMs, but with a moving threshold.
Of course, linear models do not accurately describe the neuronal data from experiments on DM. However, previous work has shown that attractor network models for twoalternative DM operate in the vicinity of pitchfork bifurcation, which is what underlies the winner-take-all competition leading to the decision dynamics [17]. In this regime the neuronal dynamics is well described by a stochastic normal-form equation which right at the bifurcation is precisely equivalent to the DDM with an additional cubic nonlinearity. This nonlinear DDM fits behavioral data extremely well, including both correct and error reaction times. In the second part of the paper I will show how such normal-form equations can be derived for an arbitrary number of neuronal populations. These equations can be thought of as nonlinear DDMs and, in fact, are identical to the linear DDMs with the addition of quadratic nonlinearities (for n > 2). Amazingly, the dynamics of such a nonlinear DDM can be recast as the diffusion of particle on a potential, which is obtained analytically, for arbitrary n.

Results
The canonical drift diffusion model (DDM) can be written where X is the decision variable, and μ is the drift or the degree of evidence in favor of one choice over the other: we can associate choice 1 with positive values of X and choice 2 with negative values. The Gaussian process ξ (t) represents noise and/or uncertainty in the integration process, with ξ (t) = 0 and ξ (t)ξ (t ) = σ 2 δ(tt

Derivation of a DDM for two-alternative DM
The DDM can be derived from a set of linear equations which model the competition between two populations of neurons, the activity of each of which encodes the accumulated evidence for the corresponding choice. The equations are τṙ 1 = -r 1 + sr 1cr I + I 1 + ξ 1 (t), where r I represents the activity of a population of inhibitory neurons. The parameter s represents the strength of excitatory self-coupling and c is the strength of the global inhibition. The characteristic time constants of excitation and inhibition are τ and τ I respectively. A choice is made for 1 (2) whenever r 1 = r th (r 2 = r th ).
It is easier to work with the equations if they are written in matrix form, which iṡ where In order to derive the DDM, I express the firing rates in terms of three orthogonal basis functions: one which represents a competition between the populations e 1 , a second for common changes in population rates e c and a third which captures changes in the activity of the inhibitory cells e I . Specifically I write where e 1 = (1, -1, 0) and e C = (1, 1, 0) and e I = (0, 0, 1). The decision variable will be X, while m C and m I stand for the common mode and inhibitory mode respectively. The dynamics of each of these modes can be isolated in turn by projecting Eq. (2) onto the appropriate eigenvector. For example, the dynamics for the decision variable X are found by projecting onto e 1 , namely e 1 ·ṙ = e 1 · (L L Lr) + e 1 · I + e 1 · ξ , ( 5 ) and similarly the dynamics for m C and m I and found by projecting onto e C and e I respectively. Doing so results in the set of equations If the self-coupling sits at the critical value s = 1, then these equations simplify to τ IṁI = -m I + gm c + I I + ξ I (t), from which it is clear that the equation for X describes a drift-diffusion process. It is formally identical to Eq. (1) with μ = I 1 -I 2 2 and ξ (t) = ξ 1 (t)-ξ 2 (t)

2
. Importantly, X is uncoupled from the common and inhibitory modes, which themselves form a coupled subsystem. For s = 1 the decision variable still decouples from the other two equations, but the process now has a leak term (or ballistic for s > 1) [18]. It has therefore been argued that obtaining a DDM from linear neuronal models requires fine tuning, a drawback which can be avoided in nonlinear models in which the linear dynamics is approximated via multistability; see e.g. [19,20]. If one ignores the noise terms, the steady state of this linear system is (X, m C , m I ) = (X 0 , M c , M I ) where One can study the stability of this solution by considering a perturbation of the form (X, m C , m I ) = (X 0 , M c , M I ) + (δX, δm C , δm I )e λt . Plugging this into Eq. (7) one finds that there is a zero eigenvalue associated with the decision variable, i.e. λ 1 = 0, whereas the eigenvalues corresponding to the subsystem comprising the common and inhibitory modes are given by which always have a negative real part. Therefore, as long as τ I is not too large, perturbations in the common mode or in the inhibition will quickly decay away. This allows one to ignore their dynamics and assume they take on their steady-state values. Finally, the bounds for the decision variable are found by noting that r 1 = X + M C and r 2 = -X + M C . Therefore, given that the neuronal threshold for a decision is defined as r th , we find that θ = ±(r th -M C ).

Derivation of a DDM for three-alternative DM
I will go over the derivation of a drift-diffusion process for three-choice DM in some detail for clarity, although conceptually it is very similar to the two-choice case. Then the derivation can be trivially extended to n-alternative DM for any n. The linear rate equations are once again given by Eq. (3), with One once again writes the firing rates in terms of orthogonal basis functions, of which there must now be four. The common and inhibitory modes are the same as before, whereas now there will be two distinct modes to describe the competition between the three populations, in contrast to just a single decision variable. Any orthogonal basis in the 2D space for competition is equally valid. However, in order to make the choice systematic, I assume that the first vector is just the one from the two-alternative case, namely (1, -1, 0, 0), from which it follows that the second must be (up to an amplitude) (1, 1, -2, 0). Then for n alternatives I will always take the first n -2 basis vectors to be those from the n -1 case. The last eigenvector must be orthogonal to these. Specifically, for n = 3, One projects Eq. (3) onto the four relevant eigenvectors, which leads to the following equations: When s = 1 then the first two equations in Eq. (11) describe a drift-diffusion process in a two-dimensional subspace, while the coupled dynamics of the common and inhibitory modes are once again strongly damped. The DDM for three-alternative DM can therefore be written Note that the dynamics of the two decision variables X 1 and X 2 are coupled through the correlation in their noise sources. The decision boundaries are set by noting that Therefore, given that the neuronal threshold for a decision is defined as r th we can set three decision boundaries: 1 Population 1 wins if X 2 = -X 1 + r th -M C , 2 Population 2 wins if X 2 = X 1 + r th -M C and 3 Population 3 wins if X 2 = -(r th -M C )/2. These three boundaries define a triangle in (X,Y)-space over which the drift-diffusion process take place.

Derivation of DDMs for n-alternative DM
The structure of the linear rate equations Eq. (3) can be trivially extended to any number of competing populations. In order to derive the corresponding DDM one need only properly define the basis functions for the firing rates, which was described above. The common and inhibitory modes always have the same structure. If the basis functions are e i and the corresponding decision variables X i for i = [1, n -1], then the firing rates are r = n-1 i=1 e i X i (t) and it is easy to show that the dynamics for the kth decision variable is given by as long as s = 1. The decision boundaries are defined by setting the firing rates equal to their threshold value for a decision, i.e.
where u = (1, 1, . . . , 1). The basis set proposed here for n-alternative DM is to take for the kth eigenvector e k = (1, 1, . . . , 1, -k, 0, . . . , 0), where the element -k appears in the (k + 1)st spot, and which is a generalization of the eigenvector basis taken earlier for two-and three-alternative DM. With this choice, the equation for the kth decision variable can be written The firing rate for the ith neuronal population can then be expressed in terms of the decision variables as which, given a fixed neuronal threshold r th , directly gives the bounds on the decision variables. Namely, the ith neuronal population wins when -(i -1)x i-1 + n-1 l=i x l + M C > r th . The n-alternative DDM reproduces the well-known Hick's law [21], which postulates that the mean reaction time (RT) increases as the logarithm of the number of alternatives, for fixed accuracy; see

DDM for multiple alternatives as the continuous-time limit of the MSPRT
In the previous section I have illustrated how to derive DDMs starting from a set of linear differential equations describing an underlying integration process of several competing Figure 1 Hick's Law for multi-alternative DDMs. If the accuracy is held fixed as the number of alternatives n is increased, then the mean reaction time increases as the logarithm of n. This is shown for both an accuracy of 0.8 (circles) and 0.6 (squares). The inset shows that the variance in the RT increases proportional to the mean RT. Parameter values are σ = 1, τ = 20 ms. The thresholds θ Pc=0.8 = 3 and θ Pc=0.6 = 0.9 for n = 2, and are increased for n > 2 to achieve the same accuracy. The initial condition is always taken to be X i = 0 for all i. The symbols are the average of 10,000 runs. The solid curves in the main panel are fits to the function RT = a + b ln (c + n). The solid lines in the inset are best-fit linear regressions data streams. An alternative computational framework would be to apply a statistical test directly to the data, without any pretensions with regards to a neuronal implementation. In fact, the development of an optimal sequential test for multiple alternatives followed closely on the heels of Wald's work for two alternatives in the 1940s [22]. Subsequent work proposed a Bayesian framework for a multiple sequential probability ratio test (MSPRT) [14][15][16]. Here I will show how such a Bayesian MSPRT is equivalent to a corresponding DDM in the continuum limit, albeit with moving thresholds.
I will follow the exposition from [16] in setting up the problem, although some notation differs for clarity. I assume that there are n possible alternatives, and that the instantaneous evidence for an alternative i is given by z i (t), which is drawn from a Gaussian distribution with mean I i and variance σ 2 . The total observed input over all alternatives n and up to a time t is written I and the hypothesis that alternative i has the highest mean H i . Then, using Bayes theorem, the probability that alternative i has the highest mean given the total observed input is Furthermore one has Pr(I) = n k=1 Pr(I|H k ) Pr(H k ). Given equal priors on the different hypotheses, Eq. (20) can be simplified to Finally, the log-likelihood of the alternative i having the largest mean up to a time t is L i (t) = ln Pr(H i |I), which is given by The choice of Gaussian distributions for the input leads to a simple form for the loglikelihoods. Specifically, the time rate-of-change of the log-likelihood for alternative i is given bẏ where where ξ i is a Gaussian white noise process with zero mean and variance σ 2 .
I can now write the L i in terms of the orthogonal basis used in the previous section, Formally projecting onto each mode in turn leads precisely to the DDMs of Eq. (18) (with τ = 1) for the decision variables, while the dynamics for the common mode iṡ For Bayes-optimal behavior, a choice i should be chosen if the log-likelihood exceeds a given threshold, namely if Note that the log-likelihood is always negative and hence does not represent a firing rate, as in the differential equations studied in the previous section. This does not pose a problem since we are simply implementing a statistical test. On the other hand, an important difference with the case studied previously is the fact that, for the neuronal models, the common mode was stable and converged to a fixed point. Therefore, the decision variable dynamics was equivalent to the original rate dynamics with a shift of threshold. Here, that is not the case. The common mode represents the normalizing effect of the log-marginal probability of the stimulus, which always changes in time. Specifically, if we assume that I l > I i for all i, namely that the mean of the distribution is greatest for alternative l, then the expected dynamics of the common mode at long times are Therefore the DDMs are only equivalent to applying the Bayes theorem if the threshold is allowed to vary in time. In this case they are, in fact, mathematically identical and thus give the same accuracy and reaction time distributions, as shown in Fig. 2.
On the other hand, an equivalent DDM with a fixed threshold has worse accuracy and shorter reaction times. The way I choose the parameters for a fair comparison is to set the fixed threshold such that the mean reaction time is identical to the Bayesian model for zero coherence. Another important difference is that the error RTs are longer than the correct RTs for the Bayesian model, see Fig. 2B, an effect which is commonly seen in experiment (and is also reproduced by the nonlinear DDMs studied in the next section of the paper) [17]. On the other hand correct and error RTs for the DDMs are always the same.
The nonlinear transfer function φ (φ I ) does not need to be specified in the derivation. The noise sources ξ i are taken to be Gaussian white noise and hence must sit outside of the transfer function; they therefore directly model fluctuations in the population firing rate rather than input fluctuations. Input fluctuations can be modeled by allowing for a non-white noise process and including it directly as an additional term in the argument of the transfer function. Note that here I assume the nonlinearity is a smooth function. This is a reasonable assumption for a noisy system such as a neuron or neuronal circuit. Nonsmooth systems, such as piecewise linear equations for DM, require a distinct analytical approach; see, e.g., [23]. The details of the derivation for two alternatives can be found in [17], but here I give a flavor for how one proceeds; the process will be similar when there are three or more choices, although the scaling of the perturbative expansion is different. One begins by ignoring the noise sources and linearizing Eq. (28) about a fixed-point value for which the competing populations have the same level of activity, and hence also I 1 = I 2 . Specifically one takes (r 1 , r 2 , r I ) = (R, R, R I ) + (δr 1 , δr 2 , δr I )e λt , where δr 1. In vector form this can be written r = R + δre λt . Plugging this ansatz into Eq. (28) and keeping only terms linear in the perturbations leads to the following system of linear equations: where and the slope of the transfer function φ is calculated at the fixed point. Note that the matrix Eq. (30) has a very similar structure to the linear operator in Eq. (3). This system of equations only has a solution if the determinant of the matrix is equal to zero; this yields the characteristic equation for the eigenvalues λ. These eigenvalues are Note that λ 1 = 0 if 1sφ = 0, while the real part of the other two eigenvalues is always negative. This indicates that there is an instability of the fixed point in which the activity of the neuronal populations is the same, and that the direction of this instability can be found by setting 1sφ = 0 and λ = 0 in Eq. (30). This yields L L L 0 · δr = 0, where the solution of which can clearly be written δr = (1, -1, 0). This is the same competition mode as found earlier for the linear system.

A brief overview of normal-form derivation
At this point it is still unclear how one can leverage this linear analysis to derive a DDM. Specifically, and unlike in the linear case, one cannot simply rotate the system to uncouple the competition dynamics from the non-competitive modes. Also, note that the steady states in a nonlinear system depend on the external inputs, whereas that is not the case in a linear system. In particular, the DDM has a drift term μ which ought to be proportional to the difference in inputs I 1 -I 2 , whereas to perform the linear stability we assumed I 1 = I 2 . Indeed, if one assumes that the inputs are different, then the fixed-point structure is completely different. The solution is to assume that the inputs are only slightly different, and formalize this by introducing the small parameter . Specifically, we write I 1 = I 0 + 2 I + 3Ī 1 , and I 2 = I 0 + 2 I + 3Ī 2 . In this expansion, I 0 is the value of the external input which places the system right at the bifurcation in the zero-coherence case. In order to describe the dynamics away from the bifurcation we also allow the external inputs to vary. Specifically, I represents the component of the change in input which is common to both populations (overall increased or decreased drive compared to the bifurcation point), whileĪ i is a change to the drive to population i alone, and hence captures changes in the coherence of the stimulus. The particular scaling of these terms with is enforced by the solvability conditions which appear at each order. That is, the mathematics dictates what these are; if one chose a more general scaling one would find that only these terms would remain.
The firing rates are then also expanded in orders of and written where r 0 are the fixed-point values, e 1 is the eigenvector corresponding to the zero eigenvalue and X is the decision variable which evolves on a slow-time scale, T = 2 t. The slowtime scale arises from the fact that there is an eigenvector with zero eigenvalue; when we change the parameter values slightly, proportional to , the growth rate of the dynamics along that eigenvector is no longer zero, but still very small, in fact proportional to 2 in this case. The method for deriving the normal-form equation, i.e. the evolution equation for X, involves expanding Eq. (28) in . At each order in there is a set of equations to be solved; at some orders, in this case first at order O( 3 ), the equations cannot be solved and a solvability condition must be satisfied, which leads to the normal-form equation.

The normal-form equation for two choices
Following the methodology described in the preceding section leads to the evolution equation for the decision variable X, where for the case of Eq. (28), η = φ /2, ξ (t) = (ξ 1 (t)ξ 2 (t))/2 and the coefficients μ and γ are see [17] for a detailed calculation. Equation (35) provides excellent fits to performance and reaction times for monkeys and human subjects; see Fig. 3 from [17]. It is important to note that the form of Eq. (35) only depends on there being a two-way competition, not on the exact form of the original system. As an example, consider another set of firing rate equations τṙ 1 = -r 1 + φ(sr 1cr 2 + I 1 ), where rather than model the inhibition explicitly, an effective inhibitory interaction between the two populations is assumed. In this case the resulting normal-form equation is still Eq. (35). In fact, performing a linear stability analysis on Eq. (37) yields a null eigenvector e 1 = (1, -1). This indicates that the instability causes one population to grow at the expense of the other, in a symmetric fashion, as before. This is the key point which leads to the normal-form equation. More specifically we see that for both systems r 1 = R + X while r 2 = R -X, which means that if we flip the sign on the decision variable X and switch the labels on the neuronal populations, the dynamics is once again the same. This reflection symmetry ensures that all terms in X will have odd powers in Eq. (35) [17]. It is broken only when the inputs to the two populations are different, i.e. by the first term on the r.h.s. in Eq. (35).
As we shall see, the stochastic normal-form equation, Eq. (35), which from now on I will refer to as a nonlinear DDM, has a very different form from the nonlinear DDMs for n > 2. The reason is, again, the reflection symmetry in the competition subspace for n = 2, which is not present for n > 2. Therefore, for n > 2 the leading-order nonlinearity is quadratic, and, in fact, a much simpler function of the original neuronal parameters.

Three-alternative forced-choice decision making
The derivation of the normal form, and the corresponding DDM for three-choice DM differs from that for two-alternative DM in several technical details; these differences continue to hold for n-alternative DM for all n ≥ 3. Therefore I will go through the derivation in some detail here and will then extend it straightforwardly to the other cases.
(38) I first ignore the noise terms and consider the linear stability of perturbations of the state in which all three populations have the same level of activity (and so I 1 = I 2 = I 3 = I 0 ), i.e. r = R + δre λt , where R = (R, R, R, R I ). This once again leads to a set of linear equations L L L(λ)δr = 0. The fourth-order characteristic equation leads to precisely the same eigenvalues as in the two-choice case, Eqs. (31) and (32), with the notable difference that the first eigenvalue has multiplicity two. This means that if 1sφ = 0 then there will be two eigenvalues identically equal to zero and two stable eigenvalues. This is the first indication that the integration process underlying the DM process for three choices will be two-dimensional. The eigenvectors for the DM process are found by solving L L L 0 δr = 0, where There are many possible solutions; a simple choice would be e 1 = (1, -1, 0, 0) and e 2 = (1, 1, -2, 0), and so e T 1 · e 2 = 0. Note that in this linear subspace any valid choice of eigenvector will have the property that the sum of all the elements will equal zero; this will be true whatever the dimensionality of the DM process and reflects the fact that all of the excitatory populations excite the inhibitory interneurons in equal measure.
To derive the normal form I once again assume that the external inputs to the three populations differ by a small amount, namely (I 1 , I 2 , I 3 ) = (I 0 , I 0 , I 0 ) + 2 (Ī 1 ,Ī 2 ,Ī 3 ), and then expand the firing rates as r = R + (e 1 X 1 (T) + e 2 X 2 (T)) + O( 2 ), where the slow time is T = t. Note that the inputs are only expanded to second order in , as opposed to third order as in the previous section. The reason is that the solvability condition leading to the normal-form equation for the 2-alternative case arises at third order. This is due to the fact that the bifurcation has a reflection symmetry, i.e. it is a pitchfork bifurcation and so only odd terms in the decision variable are allowed. The lowest-order nonlinear term is therefore the cubic one. On the other hand, for more than two alternatives there is no such reflection symmetry in the corresponding bifurcation to winner-take-all behavior. Therefore the lowest-order nonlinear term is quadratic, as in a saddle-node bifurcation.
I expand Eq. (38) in orders of . In this case a solvability condition first arises at order 2 , which also accounts for the different scaling of the slow time compared to two-choice DM. Note that there are two solvability conditions, corresponding to eliminating the projection of terms at that order onto both of the left-null eigenvectors of L L L 0 . As before, the left-null eigenvectors are identical to the right-null eigenvectors. Applying the solvability condition yields the normal-form equations The nonlinear DDM Eq. (40) provides an asymptotically correct description of the full dynamics in Eq. (38) in the vicinity of the bifurcation leading to the decision-making behavior. Figure 3(A) shows a comparison of the firing rate dynamics with the nonlinear DDM right at the bifurcation. The appropriate combinations of the two decision variables X 1 and X 2 clearly track the rate dynamics accurately, including the correct choice (here

A note on the difference between the nonlinear DDM for 2A and 3A DM
The dynamics of the nonlinear DDM for 2A, Eq. (35), depends strongly on the sign of the cubic coefficient γ . Specifically, when γ < 0 the bifurcation is supercritical, while for γ > 0 it is subcritical, indicating the existence of a region of multi-stability for I < 0. In fact, in experiment, cells in parietal cortex which exhibit ramping activity during perceptual DM tasks, also readily show delay activity in anticipation of the sacade to their response field [24]. One possible mechanism for this would be precisely this type of multi-stability. When I = 0, i.e. when the system sits squarely at the bifurcation, Eq. (35) is identical to its linear counterpart with the sole exception of the cubic term. For γ < 0 the state X = 0 is stabilized. In fact, the dynamics of the decision variable can be viewed as the motion of a particle in a potential, which for γ < 0 increases rapidly as X grows, pushing the particle back. On the other hand, for γ > 0 the potential accelerates the motion of X, pushing it off to ±∞. This is very similar to the potential for the linear DDM with absorbing boundaries. Therefore, the nonlinear DDM for two-alternatives is qualitatively similar to the linear DDM when it is subcritical, and hence when the original neuronal system is multi-stable.
On the other hand, the nonlinear DDM for three alternatives, Eq. (40), has a much simpler, quadratic nonlinearity. The consequence of this is that there are no stable fixed points and the decision variables always evolve to ±∞. Furthermore, to leading order there is no dependence on the mean input, indicating that the dynamics is dominated by the behavior right at the bifurcation. a The upshot is that Eq. (40) is as similar to the corresponding linear DDM with absorbing boundaries as possible for a nonlinear system without fine tuning. This remains true for all n > 2.
This also means that neuronal systems with inhibition-mediated winner-take-all dynamics are generically multi-stable for n > 2, although for n = 2 they need not be. This is due to the reflection symmetry present only for n = 2.
A linear stability analysis shows that the eigenvalues of perturbations of the state r = (R, R, . . . , R, R I ) are given by Eqs. (31) and (32), where the first eigenvalue has multiplicity n -1. Therefore the decision-making dynamics evolves on a manifold of dimension n -1. The linear subspace associated with this manifold is spanned by n -1 eigenvectors which are mutually orthogonal and the elements of which sum to zero. For n alternatives, we take n -1 eigenvectors of the form e k = (1, 1, . . . , 1, -k, 0, . . . , 0), for the kth eigenvector (again, the -k sits in the (k+1)-st spot). Therefore one can write Following the same procedure as in the case of three-choice DM and applying the n -1 solvability conditions at order O( 2 ), one arrives at the following normal-form equation for the kth decision variable for n alternatives: It is illuminating to write this formula out explicitly for some of the decision variables: where I have left off the noise terms for simplicity. Surprisingly, the n -1 equations for the decision variables in n-alternative DM can all be derived from a single, multivariate function: For the system of firing rate equations studied here the parameters a = φ and b = s 2 φ . Then the equation for the kth decision variable is simply The dynamics of the function f is given by where I have ignored the effect of noise. Therefore the dynamics of the decision variables can be thought of as the motion of a particle on a potential landscape, given by f . Noise sources lead to a diffusion along this landscape.

Discussion
In this paper I have illustrated how to derive drift-diffusion models starting from models of neuronal competition for n-alternative decision-making tasks. In the case of linear systems, the derivation consists of nothing more than a rotation of the dynamics onto a subspace of competition modes. This idea is not new, e.g. [13], although I have made the derivation explicit here, and have chosen as a model of departure one in which inhibition is explicitly included as a dynamical variable. It turns out that a Bayesian implementation of a multiple sequential probability ratio test is also equivalent to a DDM in the continuum limit, albeit with time-varying thresholds. For nonlinear systems, the corresponding DDM is a stochastic normal form, which is obtained here using the method of multiple-scales [25]. The nonlinear DDM was obtained earlier for the special case of two-alternative DM [17]. For four-alternative DM the nonlinear DDM was obtained with a different set of competitive basis functions [26] to describe performance and reaction time from experiments with human subjects. This led to a different set of coupled normal-form equations from those given by Eq. (42), although the resulting dynamics is, of course, the same. The advantage of the choice I have made in this paper for the basis functions, is that they are easily generalizable for any n, leading to a simple, closed-form expression for the nonlinear DDM for any arbitrary number of alternatives, Eq. (42).
An alternative approach to describing the behavior in DM tasks, is to develop a statistical or probabilistic description of evidence accumulation; see, e.g., [1,13,15,22,27]. Such an approach also often leads to a drift-diffusion process in some limit, as is the case for the Bayesian MSPRT studied here, and see also [28]. In fact, recent work has shown that an optimal policy for multiple-alternative decision making can be approximately implemented by an accumulation process with time-varying thresholds, similar to the Bayesian model studied in this manuscript [29]. From a neuroscience perspective, however, it is of interest to pin down how the dynamics of neuronal circuits might give rise to animal behavior which is well described by a drift-diffusion process. This necessitates the analysis of neuronal models at some level. What I have shown here is that the dynamics in a network of n neuronal populations which compete via a global pool of inhibitory interneurons, can in general be formally reduced to a nonlinear DDM of dimension n -1. The nonlinear DDMs differ from the linear DDMs through the presence of quadratic (or cubic for n = 2) nonlinearities which accelerate the winner-take-all competition. In practical terms this nonlinear acceleration serves the same role as the hard threshold in the linear DDMs. Therefore the two classes of DDMs have quite similar behavior.
The DDM is one of the most-used models for fitting data from two-alternative forcedchoice decision-making experiments. In fact it provides fits to decision accuracy and reaction time in a wide array of tasks, e.g. [2,3,30]. Here I have illustrated how the DDM can be extended to n alternatives straightforwardly. It remains to be seen if such DDMs will fit accuracy and reaction times as well as their two-alternative cousin, although one may refer to promising results from [12] for three alternatives. Note also that the form of the nonlinear DDMs, Eqs. (44) and (45) does not depend on the details of the original neuronal equations; this is what is meant by a normal-form equation. The only assumptions needed for the validity of the normal-form equations are that there be global, nonlinear competition between n populations. Of course, if the normal form is derived from a given neuronal model, then the parameters a and b of the nonlinear potential Eq. (44) will depend on the original neuronal parameters.
As stated earlier, the nonlinear DDMs can have dynamics quite similar to the standard, linear DDM with hard thresholds. Nonetheless, there are important qualitative differences between the two classes of models. First of all, both correct and error reaction-time distributions are identical in the linear DDMs, given unbiased initial conditions, whereas the nonlinear DDMs generically show longer error reaction times [17], also a feature of the Bayesian MSPRT. Because error reaction times in experiment indeed tend to be longer than correct ones, the linear DDM cannot be directly fit to data. Rather, variability in the drift rate across trials can be assumed in order to account for differences in error and correct reaction times; see, e.g., [30]. Secondly, nonlinear DDMs exhibit intrinsic dynamics which reflect the winner-take-all nature of neuronal models with strong recurrent connectivity. As a consequence, as the decision variables increase (or decrease) from their initial state, they undergo an acceleration which does not explicitly depend on the value of the external input. This means that the response of the system to fluctuations in the input is not the same late in a trial as it is early on. Specifically, later fluctuations will have lesser impact. Precisely this effect has been seen in the response of neurons in parietal area LIP in monkeys in two-alternative forced-choice decision-making experiments; see Fig. 10B in [31]. Given that network models of neuronal activity driving decision-making behavior lead to nonlinear DDMs, fitting such models to experimental data in principle allows one to link behavioral measures to the underlying neuronal parameters. In fact, it may be that the linear DDM has been so successful in fitting behavioral data over the years precisely because it is a close approximation to the true nonlinear DDM which arises in neuronal circuits with winner-take-all dynamics.
where to place the system at the bifurcation to winner-take-all behavior I choose 1sφ = 0. Then we have the linear operator and Given the solution r 1 , the forcing term has the form This forcing term contains a projection in the left-null eigenspace of L L L 0 . This projection cannot give rise to a solution in r 2 and so must be eliminated. The solvability conditions are therefore e 1 · (L L L 1 r 1 ) = e 1 · N 2 , e 2 · (L L L 1 r 1 ) = e 2 · N 2 , which yield Eq. (40). Note that the resulting normal form has detuning or bias terms, proportional to differences in inputs; these would be called 'drifts' in the context of a DDM, but no term linearly proportional to the normal-form amplitudes, or 'decision variables' , as in the normal form for two-choice DM. Such a term describes the growth rate of the critical mode away from the bifurcation point: attracting to one side and repelling to the other. To derive these terms we will need to continue the calculation to next order. First, we must solve for r 2 . We note that, while r 2 cannot have components in the directions e 1 or e 2 , the components in the directions e C = (1, 1, 1, 0) and e I = (0, 0, 0, 1) (the common and inhibitory modes which we discussed in the section on DDMs) can be solved for. Therefore we take r 2 = e C R 2C + e I R 2I and project e c · (L L L 0 r 2 ) = e c · N 2 , e I · (L L L 0 r 2 ) = e I · N 2 , which yields where sφ (X 1 + X 2 )((scgφ I )R 2C +Ī 1 ) sφ (-X 1 + X 2 )((scgφ I )R 2C +Ī 2 ) -2sφ X 2 ((scgφ I )R 2C +Ī 3 ) 0