{'fraction_non_alphanumeric': 0.03944343617664852, 'fraction_numerical': 0.019782214156079853, 'mean_word_length': 4.444993412384717, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dopamine firing plays a dual role in coding reward prediction errors and signaling motivation in a working memory task

Significance Recent studies have confirmed the role of dopamine firing in reward prediction error, even under perceptual uncertainty. However, little is known about dopamine behavior during the use of working memory or its role in motivation to work for reward. Here, we investigated these issues in a discrimination task. Fast dopamine responses reflected a perceptual bias while remaining consistent with the reward prediction error hypothesis. When the bias increased task difficulty, motivation positively correlated with both performance and dopamine activity. In addition, dopamine slowly ramped up in a motivation-dependent way during the working memory period. Characterizing dopamine neurons’ activity during tasks in which motivation influences behavior could importantly advance our knowledge of dopamine roles in effortful control.


Discrimination Task
This study was performed on two male monkeys, Macaca mulatta, 5-7 kg. The sensory discrimination task used here has been described previously (1, 2). The schematic representation of this task is depicted in Fig. 1A. The monkey sat on a primate chair with its head fixed. The right hand was restricted through a half-cast and kept in palm-up position. The left hand operated an immovable key (elbow at ~90°) and two push buttons in front of the animal, 25 cm away from the shoulder at eye level. The centers of the switches were located 7 and 10.5 cm to the left of the midsagittal plane. In all trials, the monkey first placed the left hand on the key, and later projected to one of the two switches. Trials began when the mechanical stimulator is lowered, indenting the fingertip of one digit of the restrained hand (Probe Down, PD). The monkey places its free hand on an immovable key (Key Down, KD). The time lag between PD and KD constitutes the reaction time (RT, Fig. 1A). After the KD, a variable delay of 1.5-3s is presented to avoid anticipatory activity before the arrival of the stimulus, followed by the first stimulus (f1), lasting 0.5s. The second stimulus (f2) is presented after a 3s delay, also lasting 0.5s. The offset of f2 signals the monkey to release the key (Key Up, KU), and report its decision by pressing one of two push buttons (PB) with the left hand (lateral push button for f2>f1, medial push button for f2<f1). Immediately after the decision report, correct discriminations were rewarded with a few drops of liquid, while incorrect discriminations received a few seconds of delay before the beginning of the next trial. Stimuli were delivered to the skin of the distal segment of one digit of the restrained right hand, via a computer-controlled stimulator (BME Systems; 2 mm round tip). Initial probe indentation was 500 μm. Vibrotactile stimuli were mechanical sinusoids pulses lasting 20ms each. Stimulation amplitudes were adjusted to produce equal subjective intensities (2). Performance was quantified through psychometric techniques (Fig. 1B, D). Animals were handled in accordance with standards of the National Institutes of Health and Society for Neuroscience. All protocols were approved by the Institutional Animal Care and Use Committee of the Instituto de Fisiología Celular (UNAM).

Recordings
Recordings were obtained with quartz-coated platinum-tungsten microelectrodes (2 to 3 MΩ; Thomas Recording) inserted through a recording chamber located over the central sulcus, parallel to the midline. Midbrain dopamine neurons were recorded in and around the substantia nigra, similar to other studies in monkeys (3,4). DA neurons were identified on the basis of their characteristic regular and low tonic firing rates (1-10 spikes per second) and by their long extracellular spike potential (2.4ms ± 0.4 SD). We furthermore verified that the 22 cells used for the study did show a positive activation to reward delivery in correct (rewarded) trials and with a pause in error (unrewarded) trials. A similar criterion has been adopted in many electrophysiological studies of midbrain DA neurons (3).

Analysis of behavioral data
Animals performed the task for multiple sessions composed of about 120 trials. Behavioral data were obtained on average from 2226 trials per stimulus class (Fig. 1D). To classify trials according to the RT we defined short-RT trials as trials with RT below the median and long-RT trials those with RT above the median (Fig. 1E-F).

Analysis of the firing rate activity
The responses (z-scores) to f1 and f2 in Fig. 2A-B were standardized with respect to a temporal window preceding the onset of the base stimulus that lasted 500ms and was centered 1000ms after KD. To estimate the temporal profile of the z-score (Fig. 4, Fig. 5A-B left and center, Fig.  5C, and Fig. 6A) we calculated the firing rate for each neuron in 250ms sliding windows shifted every 50ms and standardized it as in Fig. 2A-B. Finally, the responses (z-scores) during the delay period ( Fig. 5A-B right) were measured during its entire duration (3s) and standardized as in Fig. 2A-B.

Latency values
The firing rate time-course of the responses to f2 depended on trial uncertainty and trial outcome (Fig. 4). To determine the time of divergence between the two time-courses, we applied a receiver operating characteristic (ROC) curve analysis in each sliding temporal window. This was done within the period lasting from 300ms before f2 onset to 200ms after its offset. For each neuron we obtained the normalized firing rate (z-score) in sliding windows of 250ms shifted in 10ms steps. We used the z-scores of all neurons and trials to calculate the ROC curve at each time bin. The area under the ROC curve (AUROC) was used as the index indicating differential neuronal activity across different trial types. Values of the AUROC higher or lower than 0.5 indicated that, at the population level, one type of trial evoked a higher or lower DA response than the other. To determine the statistical significance of the computed AUROCs, we used a permutation test with 1000 resamples. Significance was determined with p<0.05 in 5 consecutive windows and the latency was defined as the first window that met this criterion. A similar analysis was used to determine significant differences in the temporal profile of the normalized activity for trials sorted according to the RT.

Dependence of the DA activity on f 1 and class number
In order to search for f1-dependent activity we performed two different tests. We first used a linear and a sigmoidal regression analysis (Fig. 6E) to assess whether the activity was monotonic with respect to f1. Then, we used a one-way analysis of variance (ANOVA) test to identify any general, non-monotonic relationship between firing rate activity and f1. We focused both analyses on the f1-stimulation period and WM delay between f1 and f2. We calculated a mean time-dependent z-score (standardized as in Fig. 2A-B) using a sliding window of 250ms moving in steps of 10ms, from 0.5s before f1 onset up to 0.5s after f2 offset. Window times with a significant monotonic signal (slope different from zero, p<0.01, for either a linear or sigmoidal fit with Q>0.05) were marked as "significantly monotonic." Window times where the ANOVA was significant (p<0.05) were marked as "significantly dependent" (1, 5). We then divided the 5s period (from 0.5s before f1 onset up to 0.5s after f2 offset) into 10 non-overlapping intervals of 500ms and counted the number of windows that were significantly linear in each interval. For each interval, we said that the f1 dependence was significantly linear if more than the 40% of windows were significantly linear, and significantly dependent if more than 40% of windows gave a significant ANOVA p-value. Fig. S3A shows the temporal evolution of the p-value resulting from these multiple ANOVA tests. Fig. S3B shows the z-scores in windows in which the dependence was significant. A similar procedure based on the ANOVA test was employed to calculate how the responses depended on the class number. The temporal evolution of the p-value resulting from the ANOVA tests and the zscores in windows with significant dependence are shown in Fig. S3C-D, respectively.

Correlations between DA and RT
To obtain the correlations between RTs and DA activity (z-scores) in Fig. 6B-C we employed data from all correct trials independently of the class presented. The z-scores were obtained as in Fig. 2A-B with specific temporal windows for each analyzed event. 250ms after cue presentation (PD), the z-score was obtained using a window of 300ms of width. Similarly, correlations at f1 presentation were computed using a window of 480ms of width centered at 240ms after the onset of the stimulus. The z-score used to obtain the correlations during the delay period (Fig. 6C left) was in a window centered 2.5 seconds after the first stimulus offset and of 1 second width. The correlation coefficient between RT and DA was calculated using the MATLAB function "corrcoef". To verify the significance of the correlations, we performed a permutation test (Fig. 6B-C left) using 10000 randomly shuffled samples. Corrcoef function also provides a p-value for the coefficient. This p-value coincided up to the second decimal point with the one obtained with the permutation test. During the 3s delay period, the temporal evolution of correlations (Fig. 6C right) was studied by obtaining the z-score in 17 windows of width 300ms evenly spaced between 0.5 and 3 seconds after the offset of f1. Correlations, which became significant towards the end of the delay, were consistently negative regardless of the number of windows as well as the width of them. Significance was assessed using corrcoef function in MATLAB.

Bayesian model for the discrimination task
The discrimination task was modeled using a Bayesian framework. The prior probabilities of f 1 and f 2 were taken to be uniform and denoted as P(f 1 i ) (i = 1, …, 6) and P(f 2 j ) (j = 1, …, 10), respectively. It was assumed that the animal knew the class structure used in the experiment (Fig. 1B-C The belief state B(k | o1 * , o2) about the class ck = (f1 k , f2 k ) (k=1, …, 12; class labels, Fig. 1C) was defined as the set of the posterior probabilities P(f1 = f1 k , f2 = f2 k | o1 * , o2) that the class ck had been presented in this trial, conditioned to the observations o2 and o1 * . It can be written as: The first factor in the above equation is a normalization factor. The second factor is a transition matrix relating the first and the second stimulation frequencies divided by the prior probability of the second stimulation frequency. Since we assumed that the animal had perfect knowledge of the class structure, the only non-zero matrix elements correspond to the 12 classes ck of the experiment. Furthermore, since f2 can only take two possible values for a given value of f1, all non-zero transition probabilities are 0.5. The last two factors then become the beliefs b1(f1 k ) and b2(f2 k ) about the first and second frequencies being those in class ck.
The sums of the B(k |o1 * ,o2)'s over classes k for each choice, gives the two belief values (f1>f2 or f1<f2). These two sums were denoted by b(H) (higher choice, f1>f2) and b(L) = 1 − b(H) (lower choice, f1<f2), and choices were made according to the larger of these two, denoted as bc. The performance can be measured by the fraction of trials of a given class in which the decision is correct.
The two unknown model parameters, the standard deviations σ1 and σ2, were fitted by minimizing the mean squared error between the model performance and the animal's performance. We found the optimal parameter values: σ1 =5.5Hz and σ2 =3.2Hz (see next section for more details about the model fitting procedure).
The statistical uncertainty U of a given trial is U=1-max[b(H), b(L)]=1-bc, which is bounded between 0 and 0.5. The uncertainty of a given class was defined as the average of the value U across trials of each class. The uncertainty in hits (or errors) of a given class was defined as the average over the correct (or wrong) trials from each class.

Bayesian model fitting procedure
Parameter fitting in Fig. 3B and S2 were made using the Simulating Annealing solver in MATLAB. For each RT condition in Fig. S2, 200 fits were performed. The two fits that yielded the lowest error were used to model the monkey's performance in Fig. 1. Given that 200 adjustments were performed, a distribution for each of the two parameters was available. The mean value for σ1 obtained in the short-RT condition was found to be significantly lower than the noise parameter for the long-RT group (one-tailed two sample t-test; p<0.001). The same result was found for the noise parameter σ2 (one-tailed two sample t-test; p<0.001).

Reinforcement Learning Model based on Belief States
To test whether the phasic responses can be attributed to dopamine reward prediction errors (RPEs) we constructed a reinforcement learning model and checked if it was able to reproduce the observed responses. Given that in the task the relevant stimuli are only partially observable (the animal is not aware of the true value of f1 and f2) we used a belief-state temporal difference (TD) model (similar to that proposed by (6)) to compute reward expectations and simulate the RPEs signaling. This was first implemented in a Bayesian module that works similarly to the Bayesian model described above and uses the same fitted values of the two noise parameters, σ1 and σ2. This module yields the belief b1(f1 i ) (i=1, …, 6) about the value of the first frequency, the belief state B(k | o1 * , o2) about the class ck = (f1 k , f2 k ) (k=1, …, 12) at the time when the second frequency is presented and the belief bc = [b(H), 1-b(H)] about which of the two frequencies is higher. Then it transmits these inference results to a TD module that selects actions and generates reward prediction errors (RPEs).
The RL model also uses a fully observable variable pm (push movement) that represents the movement towards one of the two buttons. This variable has two possible states: py ("push yes") when the decision is f1>f2 and pn ("push no") when the decision is f1<f2. Finally, the reward function, denoted by r, is a scalar function that takes two different values for correct and incorrect decisions (see Equation 8).
At the beginning of each trial the TD module calculates the value of the first stimulus as: where Q1(i) is a set of adaptable weights. The RPE at the first stimulus is .
At the second stimulus the TD module computes the value of the class: Here the QB(k) (k=1, …,12) are another set of adaptable weights. At the second stimulus the TD module also estimates the value resulting from the comparison of the two frequencies where the are a set of adaptable weights. The gi (i=1, …, G) are convenient functions that account for the contribution of the belief to this value. The functions gi were taken as the following basis functions: if c 0.1 c 0.1 and 0 otherwise. The ci's are G=11 equally spaced centroids in [0,1] and a = π /0.1.
Given the values and the RPE at the second stimulus is: .
As in the Bayesian model, we assumed that decisions were made according to the larger of the beliefs b(H) and b(L) about which of the two frequencies was the higher. The value of the response movement when action j is selected is indicated with Vrm(j) to highlight the correspondence between the action selected and the subsequent movement.
The RPE at the response movement is: and the RPE at the delivery of reward is: , where r =1 for correct discrimination and r = -0.5 otherwise.
The RPEs , and are used to update the adaptable weights Q1(i), QB(k), , the value of the movement Vrm(j).
We assumed a discount factor 1 and used the TD(λ) algorithm to update the weights at the end of each trial. The updating rule is the following: In all the equations above represents the learning rate. The value of is the temporal interval between the relevant task events, i.e. the onset of f1, the onset of f2, the response movement and the reward.   Fig. 5E). (A) Temporal evolution of the logarithm of p-value (log(p)) resulting from multiple ANOVA tests performed to check the dependence of the z-score on the f values. We calculated 1 a mean time-dependent z-score (standardized as in Fig. 2A-B) using a sliding window of 250ms shifting in steps of 10ms and performed multiple ANOVA tests sorting the z-score according to the values of f .