## Abstract

Several integrate-to-threshold models with differing temporal integration mechanisms have been proposed to describe the accumulation of sensory evidence to a prescribed level prior to motor response in perceptual decision-making tasks. An experiment and simulation studies have shown that the introduction of time-varying perturbations during integration may distinguish among some of these models. Here, we present computer simulations and mathematical proofs that provide more rigorous comparisons among one-dimensional stochastic differential equation models. Using two perturbation protocols and focusing on the resulting changes in the means and standard deviations of decision times, we show that for high signal-to-noise ratios, drift-diffusion models with constant and time-varying drift rates can be distinguished from Ornstein-Uhlenbeck processes, but not necessarily from each other. The protocols can also distinguish stable from unstable Ornstein-Uhlenbeck processes, and we show that a nonlinear integrator can be distinguished from these linear models by changes in standard deviations. The protocols can be implemented in behavioral experiments.

## 1. Introduction

Reaction time tasks have long been used to study human decision making (Luce, 1986). One paradigm requires subjects to detect or discriminate sensory signals by making a voluntary motor response. Response time distributions obtained in such perceptual-motor tasks allow inference of cognitive information processing (Posner, 1978). With the advent of methods for recording neural activity in awake behaving animals, such tasks have been adopted by neurophysiologists, especially in studies on monkeys (Schall, 2003; Gold & Shadlen, 2007), and neuronal firing rates in several brain areas have been shown to correlate with motor responses. The lateral intraparietal area, frontal eye fields, and superior colliculus (Gold & Shadlen, 2007) all exhibit activities that ramp up over time toward a fixed level before a decision is signaled (e.g., by a saccadic eye movement in the direction of the recorded response field) (Hanes & Schall, 1996; Roitman & Shadlen, 2002; Churchland, Kiani, & Shadlen, 2008). The slopes of the ramps not only correlate with task difficulty (the harder the task, the lower the slope) but also with response time (higher slopes precede faster responses). These areas may therefore provide neural substrates for integrating sensory information toward a decision criterion before a perceptual decision is made.

Various integrate-to-threshold models have been proposed to describe both response times and neurobiological mechanisms (Smith & Ratcliff, 2004), including drift-diffusion models (Ratcliff, 1978; Mazurek, Roitman, Ditterich, & Shadlen, 2003; Smith & Ratcliff, 2004; Ditterich, 2006a; Simen, Cohen, & Holmes, 2006; Ratcliff & McKoon, 2008) and attractor neural networks with mutual inhibition (Brown & Holmes, 2001; Usher & McClelland, 2001; Wang, 2002; Wong & Wang, 2006; Lo & Wang, 2006). All share a common mechanism: accumulation or integration of sensory inputs toward a prescribed threshold, the first crossing of which determines the decision and response time. These models may nonetheless be distinguished by the details of the integration process: drift-diffusion models typically accumulate evidence at a constant rate like a biased random walk, while attractor networks have unstable or stable steady states that can, respectively, accelerate or decelerate the integration process (Usher & McClelland, 2001; Wong & Wang, 2006).

Often only behavioral data are collected for human subjects, and fits based on response times and choice accuracies are sometimes unable to distinguish among competing models (Ratcliff, Zandt, & McKoon, 1999; Ratcliff, 2006). Moreover, while cellular recordings in awake behaving animals offer direct insights into the integration process, fitting of both behavioral and neural data may still not suffice, especially when models incorporate multiple features and depend on multiple parameters. For example, a drift-diffusion model with a time-dependent urgency signal (Ditterich, 2006a; Churchland et al., 2008) may be difficult to separate from a recurrent network model with strong self-excitation (Wang, 2002; Wong & Wang, 2006). Thus far, few principled attempts have been made to tease apart different integration mechanisms.

Subthreshold electrical microstimulation of neural activities in behaving animals may provide more conclusive tests (Cohen & Newsome, 2004). Ditterich, Mazurek, & Shadlen (2003) showed that microstimulation of sensory neurons affected the speed of the decision. More interestingly, in Hanks, Ditterich, and Shadlen (2006), stimulation of cells with choice targets in the recorded neuronal response fields hastened decisions, but also reduced decision speeds when choice target directions were opposite to the recorded response fields, thus providing evidence of mutual inhibition. Perturbations need not be invasive, so they can be used in human studies: Huk and Shadlen (2005), for example, employed a brief motion pulse in the background of the primary visual stimulus; although it neither determined the choice nor influenced rewards, the pulse had a significant effect on response times.

Earlier modeling efforts (Huk & Shadlen, 2005; Wong, Huk, Shadlen, & Wang, 2007; Wong & Huk, 2008) have addressed the data of Huk and Shadlen (2005), but they employed many parameters, and direct comparisons were not made between models. Here, we conduct a more rigorous study using simpler models and seek more objective comparisons among them. We approximate four neural firing rate models as linear, scalar, stochastic differential equations (SDE) and ask how their predicted response times are affected by short piecewise-constant perturbations with varying onset times and amplitudes. By comparing changes in means and standard deviations of response times, we demonstrate that the perturbations suffice to distinguish among the models. Finally, we show that a nonlinear integrator model, qualitatively similar to that of Wong and Wang (2006), behaves much like one of the linear models.

## 2. Methods

### 2.1. Reduction to One-Dimensional Linear Integrate-to-Threshold Models.

Two-alternative forced-choice decision processes can, in essence, be modeled by two populations of excitatory neurons, each endowed with self-excitatory connections and mutual inhibition via a shared inhibitory population. Each of the three populations can then be represented in a coarsely grained firing rate model by a single unit whose state describes the population-averaged activity *r _{j}*(

*t*) of the corresponding neuronal pool (Wilson & Cowan, 1972, 1973). Here we review the further reduction of a firing rate model, under suitable hypotheses, to a one-dimensional dynamical system. (For additional details, see Brown, Gao, Holmes, Bogacz, & Cohen, 2005; Bogacz, Brown, Moehlis, Holmes, & Cohen, 2006.)

Here τ_{E} and τ_{I} are the synaptic time constants for the *E* and *I* units, and *r*_{1}, *r*_{2}, and *r _{I}* are, respectively, the activity of neural units 1, 2, and

*I*. The excitatory populations 1 and 2 are selective to stimuli 1 and 2, respectively.

*F*and

_{E}*F*are their input-output functions, and the overall input to each unit, with

_{I}*i*=

*E*or

*I*, is

*I*

_{i}=

*I*

_{recurrent, i}+

*I*

_{stimulus, i}. The decision time is the first passage time from stimulus onset to the first of

*r*

_{1}or

*r*

_{2}reaching a prescribed decision threshold, which thereby signals choice 1 or 2. Since nondecision latencies (e.g., signal transduction and motor preparation) are usually assumed to be independent of stimulus strength

*I*, we model the response time as the decision time plus a constant latency. Henceforth, we use the phrases

_{stimulus}*decision time*,

*reaction time*, and

*response time*interchangeably.

*F*,

_{E}*F*may be linearized such that equation 2.1 simplifies to Here α, β, γ are the recurrent synaptic coupling strengths for self-excitation, inhibitory-to-excitatory, and excitatory-to-inhibitory connections, and

_{I}*I*

_{0,E}(

*I*

_{0,I}) is the constant background input to all the excitatory (inhibitory) cells from outside the local circuit. Unlike Wang (2002), we exclude excitatory connections between

*r*

_{1}and

*r*

_{2}and self-inhibitory connections, retaining only essential features. If τ

_{I}≪ τ

_{E}, we can further assume that the relatively fast dynamics of

*r*equilibrates rapidly, such that

_{I}*r*≈ γ(

_{I}*r*

_{1}+

*r*

_{2}) +

*I*

_{0,I}, and equation 2.2 becomes Defining a new variable

*X*≡

*r*

_{1}−

*r*

_{2}and subtracting equation 2.4 from 2.3, we obtain where

*k*≡ (α − 1)/τ

_{E}contains the excitatory coupling strength and leak, and (

*I*

_{1}−

*I*

_{2})/τ

_{E}is proportional to the difference in inputs. The background inputs

*I*

_{0,E}and

*I*

_{0,I}and coupling strengths β and γ cancel out.

*I*

_{1},

*I*

_{2}to be time varying and including additive noise, the reduced dynamics is described by a one-dimensional SDE of the form where

*b*(

*X*,

_{t}*t*) ≡

*kX*+ (

_{t}*I*

_{1}(

*t*) −

*I*

_{2}(

*t*))/τ

_{E}and σ is the standard deviation of the noise, which is assumed to be a Wiener process with increments

*dW*drawn from a normal distribution with zero mean and unit variance. Equation 2.5 provides a general description for noisy accumulator models, as illustrated schematically in Figure 1A.

_{t}### 2.2. Linear Integrate-to-Threshold Models.

We are now in a position to introduce four linear instances of equation 2.5. When *k* < 0 (e.g., self-excitation α is relatively weak), we have a classical Ornstein-Uhlenbeck (OU) process (Uhlenbeck & Ornstein, 1930; Wang & Uhlenbeck 1945). If the decision threshold does not intervene, solutions eventually decelerate and approach a stable steady-state *X*_{steady−state} = (*I*_{1}(*t*) − *I*_{2}(*t*))/(|*k*|τ_{E}). We denote this stable OU process SOU; in contrast, if *k* > 0 (stronger self-excitation), sample paths accelerate away from a fixed point in an unstable OU process, denoted UOU. SOU and UOU provide reduced descriptions of the leaky competing accumulator of Usher and McClelland (2001) that respectively produce recency and primacy effects.

*k*= 0) and constant stimuli

*I*

_{1},

*I*

_{2}, equation 2.5 becomes a pure drift-diffusion (CD) equation (Ratcliff, 1978), in which the constant drift rate

*b*=

*b*

_{0}≡ (

*I*

_{1}−

*I*

_{2})/τ

_{E}is proportional to the difference in stimuli. This CD model is similar to, but differs in detail from, the balance between leakage and inhibition in Bogacz et al. (2006).)

*b*(

*t*) =

*b*

_{0}

*t*, since similar results obtain for other functions of time, but in appendix B, we treat a modified TD model with time-varying perturbation amplitude.

### 2.3. A Nonlinear Integrate-to-Threshold Model.

*b*(

*t*) ≡

*b*

_{0}(constant) and ε represent the biased and nonbiased stimulus inputs respectively (see Figure 1B). τ

_{X}controls the overall temporal dynamics. Equation 2.8 may be derived by normal form theory (Guckenheimer & Holmes, 1983; Roxin & Ledberg, 2008) (see Figure 1B for an illustration of the branches of stable and unstable fixed points in the noise-free limit). If |ε| and |

*b*| are sufficiently small, this model without the

*X*

^{5}

_{t}term approaches those of Wang (2002), Wong and Wang (2006), and Roxin and Ledberg (2008). We include −

*X*

^{5}

_{t}to prevent unrealistic runaway activity and, for simplicity, first assume a symmetrical system with

*b*

_{0}= 0 and then consider biased stimuli |

*b*

_{0}| > 0. This model is denoted SPB.

### 2.4. Perturbation Protocols.

#### 2.4.1. Single Pulse Perturbation with Varying Onset Time.

*b*

_{1}(

*t*), under which equation 2.5 becomes

*dX*= [

_{t}*b*(

*X*,

_{t}*t*) +

*b*

_{1}(

*t*)]

*dt*+ σ

*dW*, where

_{t}*b*

_{1}(

*t*) is applied from time

*T*to

*T*+ Δ

*T*and

*b*

_{1}(

*t*) ≡ 0 otherwise. We employ two piecewise-constant forms, the first being the step function used by Huk and Shadlen (2005): in which the amplitude

*p*can be positive or negative, assisting or opposing the unperturbed drift rate

*b*

_{0}due to the primary stimulus. We fix the duration Δ

*T*at 10% of the unperturbed mean first passage time τ

_{0}, apply

*b*

_{1}at different onset times

*T*, and ask how the mean and standard deviation of the decision time change.

#### 2.4.2. Zero-Effect Pulse-Antipulse Perturbation.

*b*

_{1}leaves the mean first passage time unchanged, we call it a

*zero-effect perturbation*(ZEP). This protocol, a variant of the paired pulse suggested in Wong et al. (2007), is applied only early in the integration process to avoid interference from the decision threshold. We shall seek critical values of λ for which ZEPs occur for each model. Figures 1C and 1D illustrate both protocols.

### 2.5. Simulations and Parameter Values.

Except for the SPB model, we first consider integration to a single decision threshold in the direction of the drift. This reduces the number of parameters, simplifies mathematical analyses, and helps isolate key effects. It applies to easy tasks in which drift rates are (relatively) high and errors rare. We subsequently relax this condition in our simulations. For the SPB model, we first set *I*_{1} = *I*_{2}, corresponding to difficult tasks and high error rates, and employ two thresholds. We then consider *I*_{1} ≠ *I*_{2} and a single threshold, representing easy tasks.

Signal-to-noise ratios were chosen such that the variance in first passage (decision) times is significant, but not so great that an unreasonably large number of trials is needed to average out the noise, and we selected perturbation amplitudes, durations, and onset times such that both protocols cause small but significant effects (e.g., maximum changes of approximately 10% in mean decision times). In particular, durations were an order of magnitude smaller than mean decision times τ_{0}, as in Huk and Shadlen (2005). For the first protocol, onset times *T* were varied from stimulus onset at *t* = 0 until there were no significant effects on mean decision time. In the second protocol, we require *T* ≪ τ_{0}. We set *X*_{0} = 0 at *t* = 0 to represent unbiased initial conditions, and the remaining parameters are chosen to ensure realistic behavior (e.g., the stable fixed point is above threshold for SOU; the unstable fixed point is below *X*_{0} = 0 for UOU). Table 1 lists parameter values.

Model . | b(X, _{t}t)
. | Parameters . | ζ . | ΔT
. | p
. |
---|---|---|---|---|---|

CD | b_{0} + b_{1} | b_{0} = 5, σ = 2.449 | 20 | 0.4 | 5 |

TD | b_{0}t + b_{1} | b_{0} = 4, σ = 2.828 | 20 | 0.1 | 4 |

SOU | kX + _{t}b_{0} + b_{1} | k = −1, b_{0} = 8, σ = 1.414, | 7 | 0.4 | 2 |

UOU | kX + _{t}b_{0} + b_{1} | k = 0.2, b_{0} = 5, σ = 1.414, | 20 | 20 | 1 |

SPB | εX + _{t}X^{3}_{t} − X^{5}_{t} + b_{0} + b_{1} | ε = −0.3, 0.05, b_{0} = 0or0.004σ = 0.01, τ_{X} = 20 | ±0.75 | 5 | 0.005 or 0.0008 (b_{0} = 0.004) |

Model . | b(X, _{t}t)
. | Parameters . | ζ . | ΔT
. | p
. |
---|---|---|---|---|---|

CD | b_{0} + b_{1} | b_{0} = 5, σ = 2.449 | 20 | 0.4 | 5 |

TD | b_{0}t + b_{1} | b_{0} = 4, σ = 2.828 | 20 | 0.1 | 4 |

SOU | kX + _{t}b_{0} + b_{1} | k = −1, b_{0} = 8, σ = 1.414, | 7 | 0.4 | 2 |

UOU | kX + _{t}b_{0} + b_{1} | k = 0.2, b_{0} = 5, σ = 1.414, | 20 | 20 | 1 |

SPB | εX + _{t}X^{3}_{t} − X^{5}_{t} + b_{0} + b_{1} | ε = −0.3, 0.05, b_{0} = 0or0.004σ = 0.01, τ_{X} = 20 | ±0.75 | 5 | 0.005 or 0.0008 (b_{0} = 0.004) |

Notes: Drift-diffusion model with constant drift rate (CD) and time-varying drift rate (TD); SOU (UOU) stable (unstable) Ornstein-Uhlenbeck processes; SPB: nonlinear model. See equations 2.5 to 2.10 for details. Here ζ is the distance from starting point *X* = 0 to decision threshold, and for SPB ε = −0.3 before the stimulus appears at *t* = 0 and ε = 0.05 for *t* ⩾ 0. Perturbation durations Δ*T* and amplitudes *p* are values used in the first perturbation protocol. Parameters for the linear models with smaller signal-to-noise ratios are specified in section 3.3.

We use a forward Euler-Maruyama scheme (Higham, 2001), integrating a sample path until it hits the prescribed decision threshold and recording the corresponding first passage time, at which the trial ends and the next begins. After collecting an appropriately large ensemble, we extract the first and second moments of the first passage time from the samples. The simulation is run once for each of the different perturbations. For the models CD, TD, SOU, and UOU, the step size is 10^{−3} of a time unit and the sample size is *N* = 10^{6}, so that errors in the moments are of order . Since changes in moments due to the perturbations are of order 10^{−2}, these choices of time step and sample size reliably capture the effects of perturbation. The mean exit time for the unperturbed SPB model is about 80 time units, so that a time step of 0.01 and a sample size is 10^{6} suffices. We ran noisy simulations of this model with ε = −0.3 for *t* < 0 to represent prestimulus activity and then switched to ε = 0.05 for *t* ⩾ 0 so that *X* = 0 undergoes a pitchfork bifurcation and becomes unstable, thus forcing a choice (see Figure 1B; magnitudes of ε and σ for SPB are much smaller than the corresponding *b*_{0} and σ for the linear models because the +*X*^{3}_{t} term accelerates solutions toward the thresholds).

The parameters are not independent, and their number can be reduced by one in all models by dividing equation 2.5 by *b*_{0} and rescaling σ → σ/*b*_{0}, *p* → *p*/*b*_{0}, ζ → ζ/*b*_{0}, *k* → *k*/*b*_{0}, and the dynamical variable *X _{t}* →

*X*/

_{t}*b*

_{0}. The resulting CD and TD models are described by two parameters, while the SOU and UOU models require three and the SPB model four. Two additional parameters describe the amplitude

*p*and duration Δ

*T*of the perturbation.

We focus on qualitative patterns of changes in response times, and to make unbiased comparisons, we adopt a dimensionless measure of the changes in means and standard deviations by normalizing the relative changes with respect to their unperturbed values.

## 3. Results

### 3.1. Single Pulse Perturbations of Linear Models.

Using a single pulse perturbation with varying onset time *T*, we estimated changes in mean decision time for the four linear integrate-to-threshold models as described above, obtaining the results shown in Figure 2. Imposed early during the integration process, a brief pulse in the same (opposite) direction as the drift rate significantly decreases (increases) mean decision times in all four models, as shown by the dashed (solid) curves. These effects fade as *T* increases due to a thresholding effect: for later perturbations, more trials have already crossed threshold, and so there are relatively fewer trials that are perturbed than unperturbed. Thus, when all trials are averaged, the influence of perturbations is progressively reduced. Note that unlike the work of Huk and Shadlen (2005) and Wong et al. (2007), trials that have crossed threshold before and during perturbation are not excluded in the averaging process. This helps to reduce noise in the data, especially with late perturbation onset times.

In appendix A, we prove that this basic pattern—positive pulses advance mean decision times and negative pulses retard them—must hold provided that the perturbation occurs sufficiently early. The proof applies to a broad class of nonlinear systems including the SPB model of equation 2.8, provided that ε > 0 and thresholds lie inside the region in which drift magnitude increases with *X* (i.e., well below the activation levels of the stable fixed points).

Note that the drift-diffusion models (CD, TD) and the OU models (SOU, UOU) respond distinctly to early perturbations—the former exhibiting almost constant changes in mean decision times (see Figures 2A and 2B) and the latter respectively showing increasing and decreasing effects (see Figures 2C and 2D). Moreover, due to thresholding, SOU is alone in having an onset time for which the perturbative effect is maximized.

We next investigate how standard deviations are affected by perturbation onset time (see Figure 3). Richer patterns arise than those of Figure 2. Both CD and TD exhibit similar near-constant changes for early onsets, but the directions of changes for TD are opposite those for CD and the curves cross before reaching maxima. Standard deviations for SOU have a similar pattern to its mean, but UOU exhibits initial decreases followed by increases to a peak, unlike its mean. Optimal onset times exist for which the four standard deviations are maximally affected, and all four cases show distinct patterns.

In the TD model considered above, we assumed that the perturbation is not affected by the time-dependent gain *b*(*t*). This may not hold if the perturbation enters by the same sensory pathway as the stimuli. In appendix B, we consider a case in which perturbation and drift are affected in the same way. We show that this modified TD model differs qualitatively only for early perturbation onset, for which it yields reduced effects.

### 3.2. Zero-Effect Perturbations of Linear Models.

However, ZEPs may be analytically approximated for early onset times *T* with *T*′ ≔ *T* + Δ*T* ≪ τ_{0}, such that the probability of threshold crossing prior to *T*′ is negligible. In this case, a sufficient condition is that the PDFs of the perturbed and unperturbed processes are identical at *T*′, since for *t* > *T*′, both processes are governed by the same drift and noise. They are therefore indistinguishable in the limit *T*, Δ*T* → 0. Before deriving explicit expressions for λ* under this assumption, we present the results of numerical simulations with an intuitive explanation.

Since the unperturbed SDE (see equation 2.5) is independent of the current state *X _{t}* in the CD and TD cases, we expect that antisymmetric pulses with λ* = 1, in which the pulses precisely cancel, will produce ZEPs. In contrast, responses to inputs decay with time for SOU, implying that the first pulse must be larger than the second for their net effect to cancel at perturbation offset

*t*=

*T*′ (λ*>1). For UOU, the reverse should hold (λ* < 1). Figure 4 confirms that this is the case. Table 2 lists the parameters used in the simulations; note that the early onset time conditions are only weakly satisfied (

*T*′ ranges from 25% to 40% of τ

_{0}).

Model . | τ_{0}
. | std(τ_{0})
. | T
. | ΔT
. | p
. | λ* (thy.) . | λ* (sim.) . |
---|---|---|---|---|---|---|---|

CD | 4.005 | 0.960 | 0.5 | 0.5 | 5 | 1 | 1.0013 |

TD | 3.145 | 0.420 | 0.5 | 0.5 | 5 | 1 | 0.9978 |

SOU | 1.831 | 0.374 | 0.1 | 0.4 | 2 | 1.2214 | 1.2251 |

UOU | 2.954 | 0.142 | 0.2 | 1.0 | 2 | 0.9048 | 0.9050 |

Model . | τ_{0}
. | std(τ_{0})
. | T
. | ΔT
. | p
. | λ* (thy.) . | λ* (sim.) . |
---|---|---|---|---|---|---|---|

CD | 4.005 | 0.960 | 0.5 | 0.5 | 5 | 1 | 1.0013 |

TD | 3.145 | 0.420 | 0.5 | 0.5 | 5 | 1 | 0.9978 |

SOU | 1.831 | 0.374 | 0.1 | 0.4 | 2 | 1.2214 | 1.2251 |

UOU | 2.954 | 0.142 | 0.2 | 1.0 | 2 | 0.9048 | 0.9050 |

Note: τ_{0} and std(τ_{0}) are the first and second moments of the passage time without perturbation; *T* and Δ*T* are the onset time and total duration of the perturbation; *p* is the amplitude of the second pulse; and λ* (thy.) and λ* (sim.) are the ZEP pulse ratios predicted by the theory and obtained by averaging *p* < 0 and *p*>0 results of Figure 4.

*X*denote the unperturbed process and

_{t}*Z*the perturbed process. If

_{t}*X*

_{0}=

*Z*

_{0}at

*t*= 0, then

*X*=

_{T}*Z*at the onset time

_{T}*T*. When the perturbation ends at

*t*=

*T*′, the solution of the unperturbed SDE of equation 2.5 with constant drift rate

*b*

_{0}is while the perturbed system satisfies

*T*′], but this condition holds for

*T*′/τ

_{0}≪ 1, since the integrated drift term is almost zero at early times. Equations 3.2 and 3.3 differ in the term

*e*

^{kΔT}∫

^{ΔT}

_{0}

*e*

^{−kt}

*b*

_{1}(

*t*)

*dt*, which enters

*Z*

_{T′}due to the action of the perturbation

*b*

_{1}during the interval [

*T*,

*T*′]. The ZEP condition for arbitrary perturbations is therefore which for the piecewise constant pulses of equation 2.10 implies that

Hence, for *k* > 0 (UOU), λ* < 1, while for *k* < 0 (SOU), λ*>1. In the special case of CD, *k* = 0 and λ* = 1.

Summarizing, ZEPs occur for the CD and TD models when the opposing pulse amplitudes are equal (λ* = 1), but for SOU and UOU, λ* = *e*^{−kΔT/2} is larger and smaller than one, respectively. The ratios λ* predicted by equation 3.5 for the parameters of Table 2 are given in the final column of the table. In all four cases, they agree well with the zero crossings of the mean first passage times in Figure 4, in spite of the fact that *T*′/τ_{0} ≈ 0.25 − 0.4 is not very small.

### 3.3. Signal-to-Noise Ratio Can Influence Perturbation Effects.

To what extent do the results of section 3.1 hold when signal-to-noise ratios are reduced and a second threshold is added to track errors? To investigate this question, we select new parameter values according to the criteria of section 2.5.

For the CD model, we increase σ from 2.449 to 2.828 and set thresholds at ζ = ±5 instead of ±20, yielding an error rate of 5%. This introduces more complex behavior in which correct and error trials respond to the perturbations in different manners, as shown in Figure 5A. The combined averages (as computed in Huk & Shadlen, 2005, and Wong et al., 2007; see the black solid and dashed curves in Figure 5) preserve some features of previous results: for early perturbations, the direction of changes in mean RT agrees with Figure 2A, but the approximately constant change in mean decision time early in the trial is replaced by a monotonic decline. Increased noise yields earlier threshold crossings, advancing the thresholding effect and masking the signatures of Figure 2A. Moreover, the change in mean error RTs can reverse sign for late perturbation onsets. We checked this seemingly counterintuitive phenomenon numerically by fixing the random generating seed for noise in a sample trial, finding that perturbation of long RT trials can change an impending error to a correct choice (data not shown).

Increasing the noise level to σ = 7.071 with thresholds ζ = ±20 produces an error rate of 10% in the TD model and changes in mean RTs for correct and error choices similar to those of CD. However, averaged over correct and error choices, overall changes in mean RT exhibit a pattern similar to that of Figure 2B (see the black curves in Figure 5B). In this respect, TD maintains its signature more robustly with increasing error rate than CD, although the effects are progressively masked.

For the SOU model, we increase σ to 6.325, yielding an error rate of 17%. The optimal perturbation onset time for changes in mean RT averaged over correct and error trials is masked (compare Figure 2C with the black curves in Figure 5C). This error rate is relatively high, but we have checked that for error rates as small as 3%, SOU model features can also be masked, indicating that like the CD model, SOU is sensitive to noise.

For the UOU model, more parameters must be adjusted to obtain reasonable RT changes. Here we set σ = 2, ζ = ±10, *k* = 0.02, and *b*_{0} = 0.5, yielding a high error rate of 23.5%, but similar trends occur for error rates from 0.1% to 35%. Figure 5D shows that these are comparable to the previous UOU results with high signal-to-noise ratio. Overall, we are unable to distinguish between noise-masking effects and intrinsic features of the UOU integrator.

Using the same parameters, we find that changes in standard deviations for all the linear models show optimal perturbation onset times (see Figure 6). The CD, SOU, and UOU models are mutually indistinguishable (cf. Figures 6A, 6C, and 6D), but TD exhibits a crossing effect (see Figure 6B), similar to Figure 3B.

### 3.4. Perturbations of the Nonlinear Model.

We next apply both perturbation protocols to the nonlinear SPB model. The mean and standard deviation of decision times for the unperturbed system are 80 and 22 time units, respectively. The parameters used for the first protocol are given in Table 1, while for the ZEP, we set *p* = 0.005, *T* = 30, and Δ*T* = 5, consistent with application early in the accumulation process. Figures 7A and 7B show that the means and standard deviations of decision times behave like those for UOU (cf. Figures 2D and 3D), decreasing monotonically and approaching zero as onset times increase. ZEPs occur in the pulse ratio range 0.8 < λ* < 1 for both mean and standard deviation (see Figures 7C and 7D). Assuming a linear change in mean in Figure 7C, we find that λ* ≈ 0.87.

To compare with the behavior of the linear systems in the single threshold case, we now modify the SPB model so that it makes as few errors as they do. This is done by changing *b*_{0} from 0 to 0.004, which breaks the symmetry of the subcritical pitchfork bifurcation (see Figure 1B) so that decisions favor one choice, representing easy tasks (cf. Roxin & Ledberg, 2008); all other parameters remain the same. Figure 8 shows that the basic behavior of Figures 7C and 7D is preserved: effects on both mean and standard deviation of RTs decrease monotonically as perturbation onset time increases. While the effects on mean passage times of SPB are much like those for UOU, the fact that standard deviations of all four linear models exhibit maxima (see Figure 3) potentially allows us to distinguish SPB from them. The effects of ZEPs in this biased regime are similar to those shown in Figures 7C and 7D (data not shown).

## 4. Discussion

In this work, we have investigated how several simple stochastic integrate-to-threshold decision-making models respond to short perturbations. We focus on two-alternative forced-choice reaction time tasks for which general firing rate models with three coupled neural populations can be reduced to one-dimensional integrate-to-threshold systems. These reduced systems, which may be pure drift diffusion or Ornstein-Uhlenbeck processes, can be compared analytically and in simulations that require few parameters.

We examined two perturbation protocols, the first being a brief pulse with variable onset time and the second a pulse-antipulse pair whose amplitude ratio can be adjusted to produce minimal effects (the zero-effect perturbation, ZEP). The simulations of section 3 (see Figures 2 and 3) show that the changes in means and standard deviations of decision times allow both perturbation protocols to distinguish between drift-diffusion models and stable and unstable Ornstein-Uhlenbeck processes, provided that perturbations are delivered sufficiently early in the integration process and the signal-to-noise ratio is not too small. Changes in standard deviations can also assist in distinguishing drift-diffusion models with constant drift from those with time-dependent drift rates (cf. Figures 3A and 3B), although differences are small and may not be detectable experimentally. (As Ditterich, 2006a, notes, entire reaction time distributions can further assist in this.) For early perturbations, the ZEP conditions can be approximated analytically, making predictions that are confirmed by computer simulations (see Figure 4). While general analytical expressions for systems with finite noise appear elusive, in appendix A we prove inequalities for more general nonlinear systems that partially explain our simulation results.

Although we focus on linear models, our methods extend to other systems, and we include results for a nonlinear model that capture the dynamics near a subcritical pitchfork bifurcation, which is typical of reduced population models (Wang, 2002; Wong & Wang 2006; Roxin & Ledberg 2008). As the results of appendix A suggest, it behaves in a manner similar to the unstable OU model under both perturbation protocols, although perturbative effects on standard deviations may allow distinctions to be made (see Figures 7 and 8). We have also confirmed that the TD model's behavior remains similar to other increasing drift rates (e.g., *b*(*t*) ∼ *b*_{0}*t*^{2}; simulation results not shown).

These results reinforce the simulations and claims in Wong et al. (2007) that time-varying perturbations can in principle reveal integration mechanisms. However, some parameter ranges frustrate clear-cut comparisons among models. The first occurs for low signal-to-noise ratios and high error rates (cf. Figures 5 and 6), and so can be mitigated by using perturbations with more easily discriminated stimuli. Second, if there is a dead time between stimulus onset and the start of evidence integration, the CD and TD models can display similar signatures to the SOU model, even with high signal-to-noise ratios.

Since perturbations can be delivered through the senses as well as by direct electrical stimulation, the method can be noninvasive and is therefore appropriate for human subjects. Furthermore, unlike fitting neuronal firing rates, both perturbation protocols can distinguish linear and nonlinear UOU-type models (Wang, 2002; Wong & Wang, 2006; Roxin & Ledberg, 2008; Wong & Huk, 2008) from drift-diffusion models with an urgent time-dependent increase in drift rate (Smith et al., 2004; Ditterich, 2006a; Churchland et al., 2008). However, care is required in designing such experiments. Although a high signal-to-noise ratio can in principle help identify integrators, the reaction times can become too short, masking the signature behaviors of the integrators with the thresholding effect. A possible solution could be to prolong the response deadline in relatively simple reaction time tasks and analyze only data from long RT trials.

Perturbation inputs to decision-making circuits may originate in a single sensory pathway or in other brain areas, as in multisensory and top-down inputs (e.g., due to attention; Smith et al., 2004; Liu et al., 2008), which could modulate a decision throughout a sensorimotor pathway. The general notion of perturbation developed in the reduced model context extends to such modulatory inputs. The methods may also generalize to other tasks that recruit neural integrators (Goldman, Compte, & Wang, 2009), such as interval timing (Shea-Brown, Rinzel, Rakitin, & Malapani, 2006).

## Appendix A: Rigorous Estimates on First Passage Times

In this appendix, we analyze short single-pulse perturbations, first in the noiseless limit and then with additive noise. The inequalities proved here suggest explanations for the monotonicity of mean first passage times of SOU and UOU processes with respect to early onset times (cf. Figures 2C and 2D); moreover, they apply to more general, nonlinear one-dimensional systems such as SPB with thresholds |ζ| sufficiently close to *X _{t}* = 0.

### A.1. The Noiseless Limit.

_{[T1, T2]}denotes the indicator function that takes the value 1 for

*T*

_{1}<

*t*<

*T*

_{2}and zero otherwise). We choose and , such that

*p*and do not overlap and that neither perturbed solution reaches the threshold

*X*= ζ before . We shall compare solutions

*Y*and of the corresponding perturbed ODEs, started at the same initial condition , to show how the signs of

_{t}*f*′ and

*g*determine their threshold passage times τ and . We allow smooth nonlinear functions

*f*but assume that

*f*(

*X*) > 0 for all 0 ⩽

*X*< ζ and that

*g*does not change sign. The fact that solutions of scalar ODEs cannot cross each other is a key tool.

*Let X_{t}(x) denote the solution of equation A.1 with initial value x. Then X_{t}(x_{1}) < X_{t}(x_{2}) for any t > 0 if and only if x_{1} < x_{2}*.

This follows from the uniqueness of solutions of ODEs and the ordering of points in a one-dimensional phase space.

Since neither perturbation acts after and neither solution has reached threshold at , lemma 1 implies that if and only if . We now prove the main result summarized in Table 3, which is a corollary of proposition 1:

. | f′ ⩾ L>0
. | f′ ⩽ −L < 0
. |
---|---|---|

g>0 | ||

g < 0 |

. | f′ ⩾ L>0
. | f′ ⩽ −L < 0
. |
---|---|---|

g>0 | ||

g < 0 |

Note: *L* is some constant.

*Assume that the function g(t) in equations A.2 is strictly positive for t ∈ (0, ΔT). Then for some constant L, (a) if f′ ⩾ L>0, then for all ; (b) if f′ ⩽ −L < 0, then for all *.

*f*(

*X*)>0 and

*f*′(

*X*) ⩾

*L*>0 for

*X*⩾ 0 implies that . Letting , we have

*t*=

*T*

_{1},

*u*(

*T*

_{1}) = 0, and integrating inequality A.4 over the interval [

*T*

_{1},

*T*

_{2}], on which

*p*(

*t*) =

*g*(

*t*−

*T*

_{1}) and , yields

*p*(

*t*) ≡ 0 on , we integrate equation A.4 again to obtain where we also use inequality A.6 and the fact that the additive perturbations of equation A.2 are identical over the time intervals [

*T*

_{1},

*T*

_{2}] and to rewrite the integral term. Hence, and .

This implies that , as claimed. Proofs are similar for *g* < 0.

### A.2. Extension to Noisy Systems.

We first compare solutions *Y _{t}* and obtained from equations A.11a and A.11b under the same sample noise path with increments . Since the stochastic integral ∫

^{t}

_{0}

*dW*is almost surely continuous (Feller, 1957) and the integrated noise terms cancel, the comparison of solutions using proceeds as in the deterministic case above, and we conclude that for

_{t}*g*> 0 and

*f*′ ⩾

*L*>0, and for

*g*> 0 and

*f*′ ⩽ −

*L*< 0. Finally, averaging each process over an ensemble of sample paths, we may conclude that the mean first passage times of the two processes satisfy similar inequalities to those in the top row of table 3: and respectively. A similar argument applies to the case

*g*< 0.

## Appendix B: A TD Model with Time-Varying Perturbation Amplitude

Otherwise, we change only the perturbation amplitude, taking *p* = ±2 instead of 4. Figure 9A shows that the changes in mean and standard deviations of exit time for this model remain similar to the previous TD model, but with smaller perturbative effects near the beginning of the trial (quite similar to an SOU model) due to the monotonic time-dependence of perturbation amplitude. However, although the previous pattern of mean exit times no longer applies, the behavior remains qualitatively dissimilar to the CD, UOU, and SPB models.

## Acknowledgments

We thank C. D. Brody for suggesting a variant of the TD model, M. Usher and an anonymous reviewer for helpful suggestions that improved the letter, and Nengli Lim for pointing out the extension to noisy systems in proposition 2 of appendix A. X.Z. was supported in part by ONR grant N00014-01-1-0674, and K.F.W.L. and P.H. by PHS grants MH58480 and MH62196 and AFOSR grant FA9550-07-1-0537. The U.S. government is authorized to reproduce and distribute reprints for governmental purposes, notwithstanding any copyright notation thereon. The views and conclusions contained here are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Air Force Research Laboratory or the U.S. government.

## References

## Author notes

Note: Xiang Zhou and KongFatt Wong-Lin contributed equally to this work.