Internally generated preactivation of single neurons in human medial frontal cortex predicts volition

Abstract

Understanding how self-initiated behavior is encoded by neuronal circuits in the human brain remains elusive. We recorded the activity of 1019 neurons while twelve subjects performed self-initiated finger movement. We report progressive neuronal recruitment over ∼1500 ms before subjects report making the decision to move. We observed progressive increase or decrease in neuronal firing rate, particularly in the supplementary motor area (SMA), as the reported time of decision was approached. A population of 256 SMA neurons is sufficient to predict in single trials the impending decision to move with accuracy greater than 80% already 700 ms prior to subjects' awareness. Furthermore, we predict, with a precision of a few hundred ms, the actual time point of this voluntary decision to move. We implement a computational model whereby volition emerges once a change in internally generated firing rate of neuronal assemblies crosses a threshold.

Fig. 1

A. Schematic diagram depicting the experimental paradigm (Libet et al., 1983). Subjects were shown an analog clock and were asked to press a key with their right index finger, at will, anytime after one rotation of the clock. After the key press event (“P”), the clock dial stopped and subjects were asked to indicate the time of onset of the “urge/decision” to press the key (“W”). B–D: Distribution of W times, P times and P-W across trials and subjects. Bin size=100 ms (B, C) and 42.8 ms (D). The arrow shows the mean of the distribution (6071±3005 ms; 6264±3019 ms and 193±261 ms, mean±SD in B, C and D respectively). Medians=4964 ms, 5156 ms, 171 ms respectively. Ranges=[2795,19769] ms, [2795, 19812] ms, [43, 1455] ms respectively. W and P times are measured with respect to the trial onset time at t=0. The vertical dashed line in B and C indicates the first revolution of the clock. These distributions and mean values are very similar to those reported in earlier implementations of the same paradigm (e.g. (Haggard, 2008; Haggard and Eimer, 1999; Libet et al., 1983)). The dotted line in B and C shows an exponential fit to the behavioral data. The coarse exponential fit suggests that the response hazard function is approximately uniformly distributed (Rausand and Hoyland, 2004).

Fig. 2

A–B. Example waveforms for five single units (A) and five multi-units (B). After spike sorting, units were classified into single units or multi-units according to the criteria described in (Tankus et al., 2009). C. Distribution of the coefficient of variation of the interspike interval distribution for MUA (red) and SUA (blue). The dashed lines indicate the mean of the distribution and the horizontal bars denote one standard deviation. D. Anatomical location of electrodes in the frontal lobe displayed on a Montreal Neurological Institute (MNI) brain (average of 305 brains) (Collins et al., 1994). Each electrode included 8 recording microwires. (E–F) Raster plots and histograms showing the responses of a neuron in left ACCd displaying a significant response after W (ranksum test, p<10⁻⁶) (E), and one neuron in left pre-SMA with response onset prior to W (ranksum test, p<10⁻³) (F). All plots are aligned to W (time=0). Error bars indicate SEM (n=63 repetitions). The green line in the PSTH denotes the average time of key press across all trials. Bin size for the PSTH=100 ms.

Fig. 3

A–H. Examples of response profiles. (A–D) Neurons increasing their firing rates prior to W (p < 10⁻⁵, 10⁻⁵, 10⁻⁷, and 10⁻⁵ respectively). (E–F) Neurons decreasing their firing rates prior to W (p < 10⁻⁵, 10⁻⁴ respectively). (G–H) Neurons decreasing their firing rate prior to W and then increasing their firing rates around W (p < 10⁻³, 10⁻⁵ respectively). The conventions are as in Fig. 2E–F. (I–N). Examples of responses from several units that started to change their firing rate before the baseline period used in the text (−2500 to −1500 ms with respect to W). The responses are aligned to W (vertical black line); the vertical dashed line denotes the mean P. Only those trials where W occurred more than 5000 ms after the first turn of the clock are shown in the black trace. The dotted red trace shows all trials starting from 2500 ms before W (the black curve and the red curve do not overlap perfectly because there are more trials averaged in the red curve; the number of trials is indicated on the left of each subplot). The location of each unit is indicated in each subplot. Error bars denote SEM and are shown only every 500 ms.

Fig. 4

A. Average normalized response profile of all neurons in the frontal lobe responding prior to W, separated by whether they increase (red) or decrease (blue) their rate as W is approached (referred to as “I” and “D” respectively in the text for Increases or Decreases in firing rate). For each neuron, the baseline activity (−2500ms to −1500ms relative to W) was subtracted. For units showing increased activity before W (red), the PSTH (bin=100ms) was normalized by the maximum firing rate and for units showing decreased neural activity before W (blue), the PSTH was normalized by the minimum firing rate. Note that the responses start well before the interval used to define units as “I” or “D” (in contrast to Fig. S1G). Error bars denote SEM and are shown only every 500 ms for clarity. B. Average normalized response profile showing the temporal evolution of the responses for “I” (red) and “D” (blue) cells for MUA (left) and SUA (right). C1–C2. Average normalized firing rate of all “I” cells (C1) and “D” cells (C2) responding prior to W in each medial frontal lobe region. This plot includes both MUA and SUA (cf. B). Error bars denote SEM and are show only every 500 ms for clarity. C3/C4. Average normalized standard deviation of the firing rate of all “I” cells (C3) and “D” cells (C4) responding prior to W in each medial frontal lobe region. The format and conventions are the same as in C1–C2. For each unit, we computed the standard deviation of the firing rate across trials in each time bin and we normalized by the maximum standard deviation across all time bins. D. Percentage of frontal lobe neurons with significant change in firing rate compared with baseline (ranksum, p<0.01) as a function of time before W (Experimental Procedures). For each unit, we calculated the baseline firing rate in the window −2500ms to −1500ms relative to W (see Fig. S3E for earlier definitions of baseline period). Next, we calculated the firing rate in a 400 ms sliding window (100 ms steps) starting at time - 1500 ms to 0 ms and assessed significant changes from baseline using a ranksum test. The red and blue traces show the corresponding analyses restricted to MUA (red) and SUA (blue). The arrow indicates the percentage reported in Table 1. The horizontal dashed lines show the expected percentage (±SD) according to three different null models as described in Fig. S1 (red=“Random W”, green=“Poisson”, blue=”ISI conserved”; Fig. S1). E. Percent of neurons across brain regions with significant change in firing rate (compared with baseline) as a function of time before W.

Fig. 5

Discriminating activity from baseline on a trial-by-trial basis using a statistical classifier. A. Responses of 8 units (each in a different color) during one experimental session. Only 15 trials, randomly selected from the 53 trials in this session, are shown here for each unit. The vertical dashed line indicates the W time. B. Performance of a support vector machine (SVM) classifier in distinguishing changes in population activity with respect to baseline. At each time point t with respect to W (vertical dashed line), we considered the response of each neuron during the interval [t−200 ms;t+200 ms). We used a statistical classifier to assign the response of each neuron or each neuronal population as belonging to time t or the baseline period [−2500 ms; −2100 ms). The y-axis shows the performance of the classifier; the horizontal dashed line corresponds to chance performance obtained by random permutation of the training labels. We show the average performance level across all individual neurons in this session (gray). We next considered the entire ensemble of 37 units recorded during this experimental session (including single unit and multi-units, 22 in SMA, 8 in ACC, 7 in the medial temporal lobe). The black curve shows the performance of the classifier based on the ensemble activity; the gray shaded region indicates SEM based on 100 cross-validation steps (different random split of the data into a training set and a test set). In all cases, the reported performance levels are computed using test data not seen by the classifier during training. The two units illustrated in Fig. 2 were recorded during this session and are therefore included in the analysis.

Fig. 6

Single-trial decoding of response changes from neuronal population activity. A. Performance of the decoding classifier using a pseudo-population of varying number of units randomly sampled from the entire data set of 1019 units including both frontal and temporal regions. The horizontal dashed line indicates chance performance (50%). The red line corresponds to the classifier performance 1000 ms before W and the blue line corresponds to the classifier performance 500 ms before W as a function of the number of units used. The error bars indicate one standard deviation obtained by cross-validation from 100 random choices of the units and repetitions used for training the classifier. In all cases, the reported performance corresponds to test data not seen by the classifier during training. B. Comparison of decoding performance based on medial frontal (red) versus medial temporal (green) units (n=180 units). Note the significant advantage of medial frontal neurons over medial temporal ones. The analysis is the same as in part A except that here we select specific regions that are used to train and test the decoder. C. Comparison of the decoding performance based on 150 SMA (green), 150 pre-SMA (blue), 150 rostral ACC (red) and 150 dorsal ACC units (black). The analysis and format are the same as in part A. Note the higher classification performance of SMA over the other locations. D. Comparison of classification performance using units from the right hemisphere (green) versus units from the left hemisphere (red). The format and conventions follow the ones in part A. A population of n=268 units in each hemisphere was used (all locations combined). The horizontal dashed line shows chance performance level and the error bars were estimated by randomly shuffling the preW/baseline labels. The gray lines around the main curves show SEM over 100 cross-validation iterations. E. Comparison of classification performance using single units (red) vs. multi-units (green). A spike-sorting algorithm was used to discriminate single units (SUA) from the recorded multi-unit activity (MUA) and we automatically assigned clusters to SUA or MUA (Experimental Procedures). Here we compare the decoding performance using single-units (red, n=256) versus multi-units (green, n=256) (all locations and hemispheres combined). The format and conventions follow the ones in A. The horizontal dashed line shows chance performance level and the error bars were estimated by randomly shuffling the preW/baseline labels (one standard error over 100 cross-validation iterations). F. Comparison of classification performance using “I” cells (red, 50 units) vs. “D” cells (blue, 50 units). In this figure, all locations are combined and SUA and MUA are combined. The format and conventions are the same as in part A.

Fig. 7

Predicting the time of “urge/decision” onset (W). A–C. An SVM algorithm was used to predict the time of “urge/decision” onset (W) based on the population spiking data using 512 units. The activity of each unit was aligned to W to be able to compare activity across different recording sessions and subjects. The classifier was trained to recognize whether W had been reached or not, using windows of size 400 ms (Experimental Procedures). The binary classifier was trained using 70% of the trials and its performance was tested on the remaining 30% of the trials. The analysis window was shifted from −3500 ms up to +1000 ms with respect to W. During testing, the predicted W time was defined as the first time point when 3 out of 4 consecutive windows yielded a label indicating the occurrence of W. A. Single trial spike train example marking the position of the spikes and W. B. Spike counts in windows of size t_r=400 ms. Gray rectangles denote windows where W occurred within a time t_b ms. C. Spike count windows overlapped by 100 ms. D. Distribution of the difference between the predicted time and W (the real W corresponds to t=0 and is denoted here by the vertical dashed line). Bin size=100 ms, n=3963 trials (using cross-validation). The black and gray arrows denote the mean (−152 ms) and median (−100 ms) of the distribution respectively (standard deviation=370 ms). The dashed arrow indicates the mean value for a control case where training labels were assigned randomly (mean=1153 ms, standard deviation=995 ms). The fraction of missed trials (where the classifier could not detect W) was 2% (91% for the random label control case).