## Abstract

Self-initiated movements are reliably preceded by a gradual buildup of neuronal activity known as the readiness potential (RP). Recent evidence suggests that the RP may reflect subthreshold stochastic fluctuations in neural activity that can be modeled as a process of accumulation to bound. One element of accumulator models that has been largely overlooked in the literature is the stochastic term, which is traditionally modeled as Gaussian white noise. While there may be practical reasons for this choice, we have long known that noise in neural systems is not white – it is long-term correlated with spectral density of the form 1/*f*^{β}(with roughly 1 < β < 3) across a broad range of spatial scales. I explored the behavior of a leaky stochastic accumulator when the noise over which it accumulates is temporally autocorrelated. I also allowed for the possibility that the RP, as measured at the scalp, might reflect the input to the accumulator (i.e., its stochastic noise component) rather than its output. These two premises led to two novel predictions that I empirically confirmed on behavioral and electroencephalography data from human subjects performing a self-initiated movement task. In addition to generating these two predictions, the model also suggested biologically plausible levels of autocorrelation, consistent with the degree of autocorrelation in our empirical data and in prior reports. These results expose new perspectives for accumulator models by suggesting that the spectral properties of the stochastic input should be allowed to vary, consistent with the nature of biological neural noise.

- autocorrelation
- Bereitschaftspotential
- leaky stochastic accumulator
- readiness potential
- voluntary action

## Significance Statement

The cortical readiness potential (RP) is a gradual buildup of scalp electrical potential, and underlying neural activity in motor areas, that reliably precedes the onset of voluntary self-initiated movements by up to one second or more. More than fifty years after its discovery, the functional nature of the RP remains unclear. Here I argue, based on empirical evidence, that the RP reflects the stochastic input to an accumulation-to-bound decision process, and that this stochastic input is temporally autocorrelated, and not Gaussian white noise as it is traditionally modeled. The argument is supported by testing and confirming two novel predictions that emerge from an accumulator model when the stochastic input noise is autocorrelated rather than white.

## Introduction

Uncued voluntary movements are preceded reliably by a slow buildup of cortical activity known as the Bereitschaftspotential or readiness potential (RP; Shibasaki and Hallett, 2006). Since its discovery in the 1960s (Kornhuber and Deecke, 1965), the RP has been interpreted as a sign of movement preparation, the outcome of a preconscious neural decision to initiate an action (Libet et al., 1983). Recent evidence, however, supports a different interpretation: the RP reflects ongoing stochastic fluctuations in neural activity that favor the spontaneous emission of a movement at certain times more so than at others (Schurger et al., 2012; Murakami et al., 2014; Schmidt et al., 2016). Both the shape of the RP and the distribution of waiting times (how long the subject waits before producing a spontaneous movement) can be well described by a leaky stochastic accumulator model (Schurger et al., 2012; Murakami et al., 2014).

Accumulator models at a minimum consist of a constant term (reflecting the decision evidence) plus a stochastic term (see Materials and Methods). Integration over these “inputs” to the accumulator results in the “output” of the accumulator, commonly referred to as the decision variable because a decision is made and action initiated when the output exceeds a certain threshold. One element of accumulator models that has been largely overlooked in the literature is the stochastic term, which is traditionally modeled as Gaussian white noise. While there may be practical reasons for this assumption, it is well known that noise in neural systems is not white, but instead is pink – it is long-term correlated with spectral density of the form 1/*f*^{β}(with roughly 1 < β < 3; β= 0 for white noise) across a broad range of spatial scales (Fig. 1; Pereda et al., 1998; He et al., 2010). Whether or not accumulator models should account for the spectral properties of neural noise is an open question. At the same time, we can also ask whether or not the stochastic term in such models can be mapped onto a well-defined neural phenomenon.

A recent study found evidence for an accumulator process in area M2 of rats (homolog of human premotor cortex) performing a task where the rat could spontaneously abandon waiting for a large reward and instead opt for an immediate and certain, but smaller reward (Murakami et al., 2014). Using single-unit recordings, these authors were able to identify two different and interspersed functional subtypes of neurons: ones that behaved like the inputs to an accumulator, with tonic firing rates proportional to the length of time that the rat waited, and others that behaved like the outputs of an accumulator, with firing rates ramping up to a fixed threshold level just when the rat withdrew from the waiting station and sought out the small but certain reward. Which of these two is more likely to dominate in signals picked up by a scalp electrode, the vantage point from which the RP is commonly measured, is unknown.

Thus, here I take up two distinct but related points: (1) that accumulator models of decision-making should allow the slope (β) of the noise spectrum to vary instead of using white noise (Fig. 1), especially when the imperative (drift term) is weak relative to the noise; and (2) that a buildup in event-preceding EEG potential might reflect the average stochastic input to an accumulator process (“RP-as-input”) rather than its average output (“RP-as-output”).

With the above two points in mind, I explored the behavior of a leaky stochastic accumulator when the noise over which it accumulated was temporally autocorrelated, and empirically tested two novel predictions that emerged. I tested these two predictions on behavioral and EEG data recorded while subjects performed Libet’s self-initiated movement task (Libet et al., 1983). One prediction was that the shape of the RP should vary in a specific way as a function of the waiting time (how long the subject waits before producing a spontaneous movement). The second prediction concerned the subjective estimate of the time of the “urge” to move, what Libet et al. (1983) referred to as ‘W’ time: W time was predicted to vary as a function of waiting time, becoming earlier (with respect to movement onset) with longer waiting times. Both predictions were confirmed, with the model suggesting biologically plausible levels of autocorrelation (*β* ≅ 1.4). Taken together, the results suggest that the EEG RP may reflect the autocorrelated stochastic input to an accumulation-to-bound process, and highlight the potential utility of using biologically realistic noise in accumulation-to-bound models.

## Materials and Methods

I reanalyzed behavioral and EEG data from a prior study (Schurger et al., 2012). Details about human subjects, stimuli and task, and data preprocessing are summarized below. Details about data analyses, statistics, and computational modeling are reported in full.

### Human subjects

A total of 16 subjects participated in the study (six female, mean age 28, one left handed). Subjects were paid for their participation and all gave written informed consent. Two of the subjects did not exhibit a RP (a negative trend in electrical potential preceding movement onset at Cz or any adjacent electrode) and so were excluded from further analyses, leaving *N* = 14.

### Stimuli and task

Subjects sat ∼60 cm in front of a translucent screen onto which the stimulus was back-projected. The subject sat in a reclined position with dim ambient lighting. The stimulus was a small clock face (white on a black background) with a diameter of ∼ 6° of viewing angle. A small white dot circled smoothly around the edge of the clock dial, completing one full cycle every 3 s (equivalent to 50 ms per tick mark on the clock). Each trial began with the onset of the clock rotation. Subjects were instructed to wait for the clock to complete one full cycle. After that the subject was free to perform the instructed movement (pressing a button with the thumb of the dominant hand) at any time. Subjects were encouraged to try not to preplan the moment when they would press the button, but rather to do so spontaneously without any forethought (Libet et al., 1983). After pressing the button, the dot continued to circle the clock face for one more second and then the screen went blank. Subjects were instructed that after the screen went blank they should report (verbally) what had been the position of the dot at the moment that they first became aware of their urge to press the button. In addition to this task (the “classic” task) subjects also performed a second task which was identical except that on some trials subjects were interrupted by an auditory “click” indicating that they should press the button immediately (the “interruptus” task). The analyses reported here were performed on a combination of the data form the first task and the uninterrupted trials from the second task.

### EEG data acquisition

EEG data were recorded inside of a shielded room using a 60-channel EEG system (Elekta NeuroMag EEG/MEG system) sampled at 1000 Hz. Of the 60 channels, I used a subset encompassing the standard 10-20 montage, with the addition of electrodes C1, C2, FC1, and FC2.

### EEG data preprocessing

Data preprocessing and analysis were performed using Matlab (MathWorks Inc) along with the FieldTrip toolbox for Matlab (Oostenveld et al., 2011). Data were first downsampled to 250 Hz. Independent component analysis (ICA) was used to identify and remove ocular artifacts from the data (Jung et al., 2000), and trials with artifacts remaining after this step were excluded by visual inspection. I extracted data epochs time locked to the first button press after the start of the trial. Each data epoch covered the window of time extending from -3.5 to +1.0 s relative to the button press.

### Data analyses and statistics

#### Behavioral

I define the “waiting time” (or just “wait time”) as the amount of time, in seconds, that elapsed from the beginning of the trial until the subject first pressed the button. ‘W’ time (Libet et al., 1983) was the time, relative to the onset of the button press, that the subject reported first having been aware of the urge or decision to press the button. This was recorded by taking the clock time at which the subject reported having had the urge and subtracting it from the clock time of the actual movement. The correlation between waiting time and W time was computed using Pearson’s correlation coefficient applied both to the pooled data from all subjects and also separately for each subject. In the latter case, the resulting *r* values across subjects were tested for significance (difference from zero) using Wilcoxon's signed rank test.

### EEG

EEG data epochs were sorted according to waiting time into the lower 33rd percentile (short waiting time) and the upper 33rd percentile (long waiting time), and averaged together within each group to compare the shape of the RP for short and long waiting times. Differences in the amplitude of the RP at each time point were tested for significance using a signed rank test and then subjected to a cluster-based permutation test to correct for multiple comparisons. Unless specified otherwise, the RP was measured at electrode C1 (C2 if left handed).

#### Leaky stochastic accumulator model

Accumulator models at a minimum involve integration over an input signal (drift or “imperative”) plus Gaussian white noise. Noise in the brain, however, is temporally autocorrelated with spectral density of the form 1/*f*^{β} with typically 1 < β < 3 (“pink” noise). The model that I used was the same as that used in prior studies (Schurger et al., 2012; Murakami et al., 2014), namely the leaky stochastic accumulator model (Usher and McClelland, 2001), except that the noise was pink, with *β* allowed to vary arbitrarily, instead of white (*β* fixed at 0). Performing a run of the simulation amounts to performing numerical integration over the following differential equation:
where *I* is the drift rate (the imperative to move), *k* is leak, and ξ_{β} is noise with 1/f exponent β allowed to vary arbitrarily, from 0 (white noise) to 3; *c* is a noise scaling factor (by convention *c* = 0.1). I sometimes refer to the output of the accumulator (*x* in the equation above) as the decision variable. After each run of the simulation I extracted two “data epochs,” one from the noise input to the accumulator and one from the output of the accumulator. Both epochs were time aligned to the sample at which the output of the accumulator first crossed the threshold, and spanned the interval from 5000 samples before the crossing to 500 samples after. When I fit the event-locked input to the accumulator to the shape of the RP, I use the term RP-as-input, and when I fit the event-locked output of the accumulator to the shape of the RP, I use the term RP-as-output. The simulation is one and the same in both cases, the only difference being which by-product of the simulation (the input or the output) is considered to represent the RP.

#### Generation of pink noise

Normally, when simulating an accumulator process, a series of computer-generated Gaussian-distributed pseudo-random numbers is used to instantiate the (white) noise term. For optimal performance it is best to generate the entire time series of pseudo-random numbers all at once at the beginning of the simulation, and then step through it, rather than generating a single random number on each iteration. Having the entire time course of the noise in hand allows one to change the spectral properties of the noise as desired. Power-law or 1/*f* noise has spectral power that is inversely proportional to frequency, roughly following a negative-sloping line in log-log coordinates (He et al., 2010; Fig. 1): *P* ∝ 1/*f*^{β} with *β* commonly referred to as the “1/*f* exponent.” For the simulations reported here, I first generated a time series of Gaussian pseudo-random numbers and then altered the spectrum of this noise as follows: first the signal was converted into the frequency domain via a fast Fourier transform (*fft* in Matlab); then a given (negative) slope was imposed on the power spectrum by multiplication with a log-linear function having the desired slope; and then the signal was converted back into the time domain via an inverse Fourier transform (*ifft* in Matlab). The resulting time series, which could have any arbitrary *β* (for examples, see Fig. 1), was then used as the noise input to the simulated accumulator process. Note that I make no assumptions or claims about the spectral properties of the time series other than that they must be temporally autocorrelated. The procedure described above is simply a convenient way to obtain autocorrelated noise while parametrically varying the degree of temporal autocorrelation in the time series.

#### Data fitting

For separate runs of the simulation, the average time course of the model output or the stochastic input to the model was fit to the average RP across subjects. The time courses were time locked to the first crossing time, which is of course always determined based on the output of the accumulator. For both the RP-as-input and RP-as-output simulations the fitting procedure included a scaling factor whereby the amplitude of the simulated RP was scaled to that of the empirical RP. In each case, the normalized distribution of first-crossing times of the model was simultaneously fit to the normalized empirical waiting time distribution. Fitting was performed using multidimensional unconstrained nonlinear minimization (Nelder-Mead; function fminsearch in Matlab) on the mean-squared error between the simulated and empirical data. Appropriate starting values for the parameters, including scaling parameters, were identified by first trying a few fits by hand.

#### Modeling of ‘W’ time and derivation of the prediction

I modeled the subjective estimate of the time of the conscious urge to move (Libet’s ‘W’ time) by adding a second threshold, slightly lower than the main activation threshold (Fig. 2), similar to Ganos et al. (2015). This “advance warning” threshold is sufficiently close to the main threshold that crossing it is a good predictor that the main threshold is about to be crossed (i.e., movement is about to ensue) with very high probability. ‘W’ time in this model is simply the temporal delay between the crossing of the two thresholds, expressed with respect to the crossing of the second threshold (i.e., a negative number). Note that I do not assume that this predictive information uniquely determines ‘W’ time, but only that it informs ‘W’ time. Prior evidence indicates that neural information from both before and after movement onset can influence ‘W’ time (Sirigu et al., 2004; Lau et al., 2007; Banks and Isham, 2009; Douglas et al., 2015), and thus that the brain likely makes use of information from both before and after movement onset in its estimation of ‘W’ time. Here, I only assume that a premovement prediction about an upcoming movement and its reafferent consequences (Wolpert, 1997) is at least part of the information that contributes to the brain’s estimate of ‘W’ time.

*Estimation of the* 1/*fexponent*(β)

The 1/*f* exponent was computed as the sign-reversed slope of the best fitting line fit to the fractal component of the power spectrum in log-log coordinates. The fractal component of the power spectrum was estimated using the IRASA method of Wen and Liu (2016) and the best fitting line was estimated using ordinary least squares regression (lscov in Matlab).

#### Code accessibility

Computer code (Matlab) for running the simulations is available as Extended Data 1, and is also available online at https://bitbucket.org/aschurger/lsa_rp_model.

### Extended Data 1

Download Extended Data, ZIP file.

## Results

### Behavioral

#### Waiting time

Subjects waited on average 7.1 (±0.63 s SEM) seconds to produce a movement in the classic Libet task, and 5.4 s (±0.31 s SEM) in the interruptus task. This difference was significant (*p* < 0.001; two-tailed *t* test, *t* = 4.5, df = 14), although the two distributions were mostly overlapping (see Fig. 3*B* of Schurger et al., 2012). Since I combined noninterrupted trials from the interruptus task with data from the classic task, I performed an additional test (see below) to make sure that this average difference in waiting time did not account for any of the other results related to waiting time.

#### ‘W’ time

Libet et al. (1983) used the term ‘W’ time to refer to the subjective estimate of the time of the conscious urge to move. The mean ‘W’ time in the classic task was –142 ms (±34 ms SEM), and the mean ‘W’ time in the interruptus task was –126 ms (±44 ms SEM). The ‘W’ times in the two tasks were not significantly different across subjects (*p* = 0.45; two-sided paired samples signed rank test).

### EEG RP

Fourteen of 16 subjects exhibited the gradual negative deflection in scalp electrical potential before movement onset, at or near the vertex, that is characteristic of the RP. The two subjects that did not exhibit an RP were excluded from further analyses. The average RP across subjects (*N* = 14) is shown in Figure 3*B*. The tail of the RP extended back in time at least two full seconds before nearing its horizontal asymptote.

### Data fitting

The RP-as-output variant of the model has been shown previously to fit the data well (Schurger et al., 2012), but what about the RP-as-input variant? Although the present argument is focused on empirical predictions made by the two models, nevertheless it was important to confirm that the RP-as-input model was capable of fitting the data, before going on to test empirical predictions derived from it.

The RP-as-input model is in fact capable of producing a very good fit to the data. Figure 3 shows the normalized distribution of first crossing times of the output of the accumulator fit to the normalized empirical waiting time distribution (*n* = 14). Also shown is the average time course of the stochastic input to the accumulator (time locked to threshold crossings in the output), fit to the average RP. The parameters for the best fit were obtained after an exhaustive search of the parameter space, and were as follows: β = 1.4, I = 0.1, k = 0.6, and threshold = 0.1256. Thus, a single set of parameters exists such that the RP is well fit by the event-locked stochastic input to the accumulator and the waiting time distribution is well fit by the first-crossing time distribution of the output. It is noteworthy that the value for β was 1.4, which is approximately the value for β found in our EEG data (1.34 ± 0.14 SEM; Fig. 4) and that has been observed in the past for EEG data (Pereda et al., 1998; reported β = 1.5).

One might ask whether or not it is necessary to use 1/*f* noise as input to the accumulator. Perhaps all that is necessary is to suppress the higher frequencies, in which case a low-pass filter applied to white noise might suffice (Fig. 1*F*). One problem with this reasoning is that this procedure could be seen as simply a very crude way of approximating 1/*f* noise, which would bring us back to the original question. Nevertheless, it is unclear a-priori how well the model would fit the data using low-pass filtered white noise instead of 1/*f* noise as input to the accumulator. I tested this possibility using a first-order Butterworth low-pass filter, varying the cutoff frequency. I used a low-order filter because it has a shallow cutoff, thereby approximating more closely the relationship between frequency and power found in 1/*f* noise (with higher-order filters the fit, in particular to the RP, only became worse). I used the same exhaustive parameter search here as I had used when running the model with 1/*f* noise.

Using low-pass filtered noise as input to the accumulator, under the RP-as-input assumption, the model could fit either the RP or the wait-time distribution, at the expense of the fit to the other, but had trouble fitting both with the same set of parameters. Notably, the best fit overall using low-pass filtered noise was significantly poorer than the best fit achieved with simulated 1/*f* noise (MSE of best fit with low-pass filtered noise 0.0271 ± 0.00066 SD, vs 0.0146 ± 0.00054 SD using 1/*f* noise; *p* < 10^{−9}, two-sided Wilcoxon rank sum test). Again, this might simply reflect the fact that a low-pass filter applied to white noise yields a relatively poor approximation of 1/*f* noise. However, as stated previously, I make no claims about the necessity of “true” 1/*f* (power-law) noise, only that the input noise should be temporally autocorrelated, although it would appear that simulated 1/*f* noise does lead to a significantly better fit overall.

### Prediction 1: relationship between the shape of the RP and waiting time

The predicted relationship between waiting time and the shape of the RP under the RP-as-input and RP-as-output assumptions is shown in Figure 5, for the parameters reported above (β = 1.4, I = 0.1, k = 0.6, and threshold = 0.1256). Although these specific parameters were the ones that resulted in the best fit under the RP-as-input assumption, the predicted relationship between RPs for short versus long wait times remained qualitatively the same regardless of the specific parameters used, as long as β was more than ∼0.5. As illustrated in the figure, under the RP-as-input assumption the model predicts that the early RP will have a lower amplitude for long versus short waiting times (Fig. 5*A*). On the other hand, under the RP-as-output assumption the model predicts the opposite: a higher amplitude early RP for long versus short waiting times (Fig. 5*B*). The empirical data clearly support the RP-as-input interpretation (Fig. 6; *p* < 0.01 for the mean amplitude over the range -1.5 to -0.5 s, two-tailed signed rank test).

To understand the relationship (Fig. 5*C*,*D*), consider that the output of the accumulator will tend to have a steeper slope on trials with a short waiting time. This translates into a lower initial amplitude, when time aligned to the first crossing time. Whereas trials with a longer waiting time will tend to have a more gradual, shallow slope, and thus a higher amplitude when time aligned to the first crossing time. The input to the accumulator, on the other hand, will tend to have a higher amplitude throughout the epoch for trials with a short waiting time (the higher its amplitude the sooner its integral will cross a given threshold value). This relationship is explained schematically in Figure 5*C*,*D*, which illustrates the relationship for the simplest case where the input is constant and the output is linear. Note that this overly simplified schematic is only intended to describe the relationship between the input and output but not their shape.

### Prediction 2: correlation between waiting time and W time

As mentioned in Materials and Methods, I assume that the temporal delay between the crossing of the advance warning threshold and the main activation threshold contributes to the subjective estimate of ‘W’ time (Fig. 2). Thus, if the trajectory of the decision variable is relatively steep, then that temporal delay will be relatively short, and so ‘W’ time will be later in time (closer in time to the onset of movement). On the other hand, if the trajectory of the decision variable happens to be relatively gradual, then the temporal delay will be longer and ‘W’ time will be earlier (farther back in time from the onset of movement). Therefore, the prediction is that waiting time should be negatively correlated with ‘W’ time, with longer waiting times predicting earlier ‘W’ times and vice versa (Fig. 7). Remarkably, this prediction is confirmed, both in the pooled data from all subjects (Fig. 8; Pearson’s correlation coefficient: *r* = -0.116, *p* = 0.0003; slope of regression line: b = -0.012), and in the aggregate of individual correlations across subjects (Fig. 9; *p* = 0.0009, two-sided Wilcoxon signed rank test on the *r*’s from each individual subject, *n* = 14; *r* was negative for 13 of 14 subjects, *p* = 0.0009 binomial).

The parameters used to generate the above prediction were the same as those used for fitting the data and predicting the shape of the RP: *β* = 1.4, I = 0.1, k = 0.6, and threshold = 0.1256. The prediction was robust to small changes in the parameters. Most importantly the prediction only held for *β* more than ∼1.2, i.e., this prediction depends on the stochastic input being pink. Thus, allowing β to vary revealed a novel prediction that would not have been revealed otherwise.

The advance warning threshold may be thought of as analogous to the processing of an “efference copy” by an internal forward model (Wolpert, 1997). When the prediction threshold is crossed, information is generated indicating that a specific movement is just about to happen with high probability. After the movement is completed, this information informs the subsequent subjectively-estimated ‘W’ time. Note again that I do not assume that the predictive signal uniquely determines ‘W’ time, but only that it informs ‘W’ time. There is strong evidence that neural information from both before and after movement onset can influence ‘W’ time, and thus that the brain likely makes use of information from both before and after movement onset in making this particular judgment (Lau et al., 2007; Banks and Isham, 2009; Douglas et al., 2015). Here, I assume that a premovement prediction about an upcoming movement (and its reafferent consequences) is at least part of the information that contributes to the brain’s *post hoc* estimate of ‘W’ time.

As mentioned before, to guard against the possibility that this correlation might be driven by the difference in mean waiting time for the classic and interruptus tasks (recall that noninterrupted trials from the interruptus task were combined with data from the classic task), I tested the pooled correlation separately for the data from the classic task and the data from the interruptus task. The correlation was negative and highly significant in both data sets (classic task: *r* = -0.1159, *p* < 0.001; interruptus task: *r* = -0.1156, *p* < 0.001), thus arguing against this possible interpretation.

## Discussion

Here, I explored the behavior of an accumulator model when the power spectrum of the input noise was allowed to vary from white (the standard for such models) to pink, i.e., allowing for varying levels of temporal autocorrelation in the noise. The 1/f exponent (β) of the noise spectrum thus becomes a parameter of the model, which can take values in the range 0 ≤ β < 3. Previous work has looked at stochastic differential equations that can generate autocorrelated noise in their output with a variable 1/*f* exponent (Ruseckas and Kaulakys, 2010), but this is different from the present work which used autocorrelated noise in the input to the process (in place of the Gaussian white noise used by convention).

With the 1/*f* exponent (β) of the input noise as a parameter, I tested the possibility that the canonical RP reflects the average event-locked stochastic input to an accumulator rather than its average output. I first showed that the model, under the RP-as-input assumption, was able to fit the data and then tested a prediction regarding the shape of the RP. Finally, I also tested a prediction regarding ‘W’ time that depends on the input noise being autocorrelated.

Intuitively one might expect that, from the perspective of a scalp electrode, the output of a neural accumulator might be too faint to be detected, because it might be computed by a small population of neurons (compare O'Connell et al., 2012). At the same time, the stochastic input to the decision process should be at least as readily measured at the scalp if it is shared among a large number of neurons spanning a functional network (Mantini et al., 2007). Whatever the case may be, the model provided a way to distinguish between these two possibilities by generating opposing predictions depending on the substrate (input or output) of the observed average (Fig. 5).

The model also incorporated a second “warning” threshold, slightly lower than the primary activation threshold, such that when the lower threshold is crossed the probability of soon crossing the primary activation threshold is high (Fig. 7). This lower threshold allows the model to account for the estimated time of the subjective urge to move (‘W’ time) in terms of the delay between the crossing of the two thresholds. This is clearly an oversimplification of the true state of affairs: prior data have implicated factors from both before and after movement onset in the subjective estimation of ‘W’ time (Sirigu et al., 2004; Lau et al., 2007; Banks and Isham, 2009; Desmurget et al., 2009; Douglas et al., 2015). However, the model is only intended to account for one such factor, originating before movement, plausibly the generation of a forward model (Wolpert, 1997). There has been a great deal of controversy surrounding the meaning of ‘W’ time and its relationship to the onset of intention (for review, see Maoz et al., 2015).

Note that the argument here is not that pink noise is superior to white noise in stochastic accumulator models. Rather it is an argument in favor of allowing β to vary (as it does in real neural systems) in accumulator models of decision-making. The resulting value for β that is suggested by the model could in principle be near zero (i.e., white noise) or it could be roughly between 1 and 2 (pink noise). The point is that having β as a parameter adds a new dimension to the model allowing it to account more fully for the neural data, by fitting its spectral properties. One advantage of doing so is that the stochastic input may then behave in a way that is amenable to modeling, e.g. the the RP-as-input variant of the model can only fit the data when β is non-zero. Another advantage is that with variable β, we can account for properties of the data that were previously overlooked, potentially leading to new predictions.

Allowing the spectral properties of the input noise to vary here led to two novel predictions. The first prediction was that the shape of the RP should vary as a function of the time elapsed between the beginning of the trial and the onset of the self-initiated movement (the waiting time). The relationship between the shape of the RP and waiting time was reversed depending on whether the RP was modeled as the average input to or average output from the accumulator (Fig. 5) thus constituting a strong test with which to adjudicate between these two possibilities. The second prediction was that the estimated time of the subjective urge to move, with respect to movement onset, should be negatively correlated with the waiting time, being earlier in time (with respect to movement onset) for longer waiting times, and vice versa. Both of these predictions were confirmed, lending support to the hypotheses encapsulated in the model. The results are important because they drive a wedge between two different possible (and plausible) interpretations of the RP: (1) as the average event-locked stochastic input to an integration-to-bound process, or (2) as the average event-locked output from an integration-to-bound process. Both results point to the former interpretation.

Note that Caspar and Cleeremans (2015) found that the amplitude of the RP was lower among subjects with late ‘W’ times (“short W group”) compared to subjects with early ‘W’ times (“long W group”), whereas Haggard and Eimer (1999) found the opposite. These studies did not group their data according to waiting time as I have done here, but one can infer the following: If longer waiting times are associated with both earlier ‘W’ times and a lower-amplitude early RP, then early ‘W’ times may be associated with a lower-amplitude early RP. This is broadly consistent with Caspar and Cleeremans (2015), although the difference here was in the early RP, whereas these two prior studies report a difference in the amplitude of the RP overall. Note also that the results of both Haggard and Eimer (1999) and Caspar and Cleeremans (2015) depend on the choice of baseline which has recently been a topic of debate (Khalighinejad et al., 2018). I did not apply any baseline correction to the EEG data reported here.

The hypotheses tested in the present study are based on two premises. One is that, when movements are made spontaneously and voluntarily, a subset of neurons in premotor areas behave like inputs to an accumulator while other neurons behave like outputs (Murakami et al., 2014). This first premise raises the question of which subpopulation of neurons is responsible for the RP, as captured from the perspective of a scalp electrode, or other mass-action recording modality. The second premise is that real biological noise in the brain is not white, but rather has (approximately) a 1/*f* spectral profile. Although we have known about the RP for a long time, attempts to model it mechanistically have come about only recently (Schurger et al., 2012). These two premises expand the horizon of stochastic accumulator models, especially in the low-SNR regime, and add a new dimension to the kinds of hypotheses that can be derived from such models.

The model developed here accounts for both the shape of the RP and W time in terms of a third variable, the time spent waiting to produce a movement. To date, this variable has been given very little attention in the literature, but was predicted to be a factor by the model. To the extent that the waiting time in a self-initiated movement task can be considered analogous to the reaction time in a perceptual decision-making task, we can then consider some of the parallels between them. From the point of view of the accumulator model, both reaction time and waiting time reflect the very same thing: the time it takes for the integral of evidence plus noise to reach a decision boundary. Both produce a gamma distribution, and both exhibit a positive linear relationship between the standard deviation and the mean (Schurger et al., 2012).

The simplifying assumption of Gaussian white noise in accumulator models has held strong for decades, and such models have proven themselves very capable of accounting for behavioral and neural data (Ratcliff, 1978; Ratcliff et al., 1999; Gold and Shadlen, 2001; Usher and McClelland, 2001; Ratcliff et al., 2004; Gold and Shadlen, 2007; Ratcliff and Starns, 2009, 2013). This begs the question of why we should further complicate things by introducing varying levels of autocorrelation in the input noise. The answer is simple: stochastic fluctuations in neural systems tend to be temporally autocorrelated (Destexhe et al., 2003; Rudolph and Destexhe, 2003; Henrie and Shapley, 2005; Garcia-Perez et al., 2007; Milstein et al., 2009; Van De Ville et al., 2010; He, 2011). A real neural integrator, given that it must integrate over fluctuations originating elsewhere in the brain, will likely integrate over autocorrelated fluctuations, and if so then the level of autocorrelation will likely vary between different cognitive states and between different individuals (He et al., 2010). Thus, the level of autocorrelation (1/*f* exponent) in neural time series is a meaningful and functionally relevant physiological variable that can be accounted for in models of decision-making by allowing the spectral properties of the stochastic input to vary.

While the use of simulated 1/f noise as input to the accumulator yields a very good fit to the data (wait-time distribution and RP) there might be other ways to accomplish this. I was able to achieve a moderately good fit to the data using low-pass filtered white noise, although the best fit was still significantly poorer than that achieved using simulated 1/*f* noise. Also, the (log-log) slope of the resulting spectrum was far from that observed in the empirical data (∼3.0 for low-pass filtered noise vs ∼1.4 for simulated 1/*f* noise). So, while it might be possible to fit the wait-time distribution and/or RP using low-pass filtered white noise as input, the spectrum of the resulting input time series is a very poor match to that observed empirically. Thus, one key advantage of using simulated 1/*f* noise versus low-pass filtered white noise is that it can account not only for the behavior and event-related potential, but also for the spectral properties of the EEG data, offering a more complete representation of the data. Still, no claim is made here as to whether or not true power-law noise is required. What is required, at a minimum, is that the stochastic input time series are temporally autocorrelated and that the degree of autocorrelation can be varied parametrically. The technique used here to generate simulated 1/*f* noise (see Materials and Methods) offers a means of doing that.

Accumulator models are commonly viewed as operating at a higher level of abstraction than neural models. Nevertheless, aspects of such models are commonly mapped onto neural phenomena (Gold and Shadlen, 2007). The output of the accumulator (*x*, the decision variable) is commonly taken to represent the firing rate of neurons involved in decision-making, and the constant input (*I*) is taken to represent the sensory evidence in the form of firing rates of sensory perceptual neurons that synapse, directly or indirectly, onto the decision neurons. The noise term on the other hand is not typically considered to map directly onto a specific neural phenomenon, but rather accounts for variability within and across trials. Here, I offer a more explicit treatment of what the noise term reflects: it reflects stochastic variability in neural activity originating from elsewhere in the brain, and as such perhaps should be modeled as having the same power spectrum as real neural noise. This spectrum is known to be pink rather than white.

One might argue that, although the use of autocorrelated noise in the model is more biologically realistic, it may add relatively little in terms of accounting for and helping to explain the data. As mentioned previously, standard accumulator models, with Gaussian white input noise, have been widely used in the study of perceptual decision-making and reaction time paradigms and have proven very effective in accounting for neural and behavioral data. While this may be true, the fact that the predictions made here emerge from the model only when the input noise is autocorrelated lends credence to the argument that the spectral properties of the noise used in such models does matter. Further research will be needed to bear this out, but intuitively it seems that it will likely matter most when the imperative (the drift term) is weak relative to the noise, i.e., in the context of decision-making under uncertainty.

## Acknowledgments

Acknowledgements: I thank Jacobo Sitt and Philippe Ciuciu for contributing Matlab code.

## Footnotes

The authors declare no competing financial interests.

A.S. was supported by the European Research Council Starting Grant 640626.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.