A hardwired neural circuit for temporal difference learning

Authors: Malcolm G. Campbell, Yongsoo Ra, Zhiqin Chen, Shudi Xu, Mark Burrell, Sara Matias, Mitsuko Watabe-Uchida, Naoshige Uchida

Abstract

The neurotransmitter dopamine plays a major role in learning by acting as a teaching signal to update the brain’s predictions about rewards. A leading theory proposes that this process is analogous to a reinforcement learning algorithm called temporal difference (TD) learning, and that dopamine acts as the error term within the TD algorithm (TD error). Although many studies have demonstrated similarities between dopamine activity and TD errors^1–5^, the mechanistic basis for dopaminergic TD learning remains unknown. Here, we combined large-scale neural recordings with patterned optogenetic stimulation to examine whether and how the key steps in TD learning are accomplished by the circuitry connecting dopamine neurons and their targets. Replacing natural rewards with optogenetic stimulation of dopamine axons in the nucleus accumbens (NAc) in a classical conditioning task gradually generated TD error-like activity patterns in dopamine neurons by specifically modifying the task-related activity of NAc neurons expressing the D1 dopamine receptor (D1 neurons). In turn, patterned optogenetic stimulation of NAc D1 neurons in naïve animals drove dopamine neuron spiking according to the TD error of the stimulation pattern, indicating that TD computations are hardwired into this circuit. The transformation from D1 neurons to dopamine neurons could be described by a biphasic linear filter, with a rapid positive and delayed negative phase, that effectively computes a temporal difference. This finding suggests that the time horizon over which the TD algorithm operates—the temporal discount factor—is set by the balance of the positive and negative components of the linear filter, pointing to a circuit-level mechanism for temporal discounting. These results provide a new conceptual framework for understanding how the computations and parameters governing animal learning arise from neurobiological components.

Introduction

The brain learns by comparing internal expectations against real outcomes and updating expectations when they are wrong. This idea—that a major function of the brain is to learn accurate predictions, and that it does so in part by computing and minimizing prediction errors—is an influential organizing principle in neuroscience^6^.

A prime example of this is the phasic activity of midbrain dopamine neurons, which canonically reflects the difference between predicted and actual reward—i.e. reward prediction error (RPE)—and acts as a teaching signal to update reward predictions (“value”)^1,3^. More specifically, phasic dopamine reflects a specific form of RPE called temporal difference RPE (TD RPE, or TD error), the error term in a reinforcement learning (RL) algorithm called temporal difference (TD) learning^7,8^ (Fig. 1a). TD learning is a cornerstone of most real-world successes of RL^8–11^.

The TD algorithm relies on a core derivative-like computation, from which it derives its the comparison of value estimates at adjacent time steps (Fig. 1a). The intuition is If value suddenly increases, then to the extent that this was predictable, previous states were undervalued, and vice versa. This empowers the algorithm to leverage its internal estimates to drive learning in the absence of reward feedback, using a temporally local computation to bridge events that are far apart in time.

A profound implication of this is that the time horizon of the TD learning algorithm—its temporal discount factor, γ—depends on the relative weighting of current and recent time in the local comparison of the value function at adjacent time steps. This suggests that microcircuit mechanisms underlying the derivative computation in the dopamine system may quantitively explain animals’ preference for immediate versus future rewards. Time discounting is a fundamental aspect of learning and decision making^12,13^ and may be a partly hardwired personality trait^14^, but its biological origins remain mysterious. The possibility of quantitatively linking dopamine circuits to time preference is tantalizing but unrealized.

While debates continue^15–18^, the canonical evidence supporting the TD RPE interpretation of dopamine is that dopamine neurons in trained animals (1) respond strongly to unexpected rewards, but less to expected rewards; (2) respond to reward-predicting sensory cues; and (3) become suppressed when expected rewards are omitted^1^ (Fig. 1b). In addition to these canonical TD RPE-like activity patterns, it was recently discovered that ramping dopamine activity, which can occur before reward delivery^19,20^, reflects the approximate derivative of a hidden function, and thus can also be interpreted as a TD error^5,21,22^ (Fig. 1b). Yet despite these and other pieces of correlational^2,4^ and behavioral^3,23,24^ evidence, the biological mechanisms that produce TD error-like activity in dopamine neurons, including the core derivative-like computation, remain unknown.

We hypothesized that TD RPE-like activity in dopamine neurons arises through bidirectional interactions between midbrain dopamine neurons and their primary output region, the striatum—a longstanding idea^25–28^ that remains unproven. Neural activity in the striatum correlates with value^29–32^ and influences dopamine signaling both directly and indirectly, via its two primary cell types, D1 and D2 medium spiny neurons (MSNs) (Fig. 1c). However, the idea that dopamine-striatum interactions generate TD RPE signals in dopamine neurons has been difficult to test experimentally, for two (1) tasks trained with natural rewards engage learning processes throughout the entire brain, not just the dopamine system, so signaling patterns in dopamine neurons in natural tasks cannot be assumed to arise specifically from dopamine-striatum interactions; and (2) dopamine neurons are themselves diverse, with distinct properties and functions depending in part on projection target^33–37^, so dopamine neurons may not all participate in value learning.

To unravel this complexity step-by-step, we started with an artificial opto-conditioning task with just two sensory stimuli (odors) and optogenetic stimulation of dopamine axons in the lateral nucleus accumbens (lNAc), a subregion of the striatum with concentrated TD RPE-like dopamine signaling^38,39^ and value-related neuronal activity^31,32^ in natural conditioning tasks. The results revealed that dopamine neurons are hardwired to automatically compute a temporal difference of input from D1-MSNs, resulting in patterns of TD error in dopamine signaling, and suggest that the balance of excitation and inhibition in this pathway sets the temporal discount factor of the TD algorithm.

Results

Opto-conditioning in lNAc generates TD error signals in dopamine activity

Mice were prepared for opto-conditioning by virally expressing red excitatory opsin ChrimsonR^40^ unilaterally in ventral tegmental area (VTA) dopamine neurons and green dopamine sensor GRABDA3m^41^ in lNAc (Fig. 1d). An optic fiber was implanted targeting lNAc in the same hemisphere, through which dopamine axons could be stimulated while simultaneously recording dopamine release using photometry without optical crosstalk (Extended Data Fig. 1a,b) or movement artifacts (Extended Data Fig. 1c). To ensure physiological relevance, optogenetic stimulation was calibrated to a natural reward (8μL droplet of water) (Fig. 1e, Extended Data Fig. 1b)^42^.

Mice underwent conditioning in which odors (CS+/CS−) preceded opto-stimulation (75% of CS+ trials) or no outcome (100% of CS− trials, 25% of CS+ trials) (Fig. 1f, Extended Data Fig. 1d,e). In contrast to a previous report^16^, lNAc dopamine in ChrimsonR mice (n=5), but not control tdTomato mice (n=3), gradually developed classic signatures of TD error (Fig. 1g–j, Extended Data Fig. 1f–i): a positive response to CS+ (“cue response”), preceding the time of predicted stimulation (Fig. 1g–i); a negative response aligned to omission of predicted stimulation (“omission dip”) (Fig. 1g–i), which reflected the timing of stimulation used during training (Extended Data Fig. 1j–l); and an expectation-dependent reduction in the opto-stimulation response (Fig. 1j).

To investigate these changes at the level of single dopamine neurons while preserving projection target specificity, we recorded spiking activity of antidromically optotagged lNAc-projecting VTA dopamine neurons in the same mice on a subset of sessions using Neuropixels probes (Fig. 1k–q, Extended Data Fig. 2). Projection target was identified by stimulating dopamine axons in lNAc while recording spikes at the cell body in VTA (Fig. 1k–m). As in the photometry data, cue responses and omission dips could be detected in the spiking activity of individual dopamine neurons (Fig. 1k–q, Extended Data Fig. 2). While statistically significant both within and across mice, the cue response and omission dip were small in both spiking and photometry signals (~5% of response to uncued water delivery, Fig. 1i,j, Extended Data Fig. 2)^16^, likely reflecting the minimal nature of our calibrated, local, unilateral dopamine axon stimulation, and/or possibly a missing factor that would normally accompany a natural reward.

Opto-conditioning was sufficient to alter Mice exhibited greater anticipatory facial movements following the CS+ than the CS− (Extended Data Fig. 3). However, behavioral and neural changes were dissociable, as neural changes were observed even in mice that exhibited no detectable change in movement following odor onset.

Opto-conditioning in lNAc generates a value-like signal in D1- but not D2-MSNs

To understand how opto-conditioning in lNAc generated TD error-like signals in dopamine neurons, we used Neuropixels probes to record extracellular spiking activity in the striatum, on a subset of late-training days in the same mice (Fig. 2a). We selected neurons with a positive response to either the CS+ or CS− (Methods) and plotted their activity on omission trials. In contrast to the biphasic response observed in dopamine neuron spiking (Fig. 1o,p) and dopamine release (Fig. 1g,h), striatal spiking was positively enhanced to the CS+ relative to the CS− throughout the trial (Fig. 2b, Extended Data Fig. 4), reminiscent of a “value”-like signal as opposed to the TD error-like patterns observed in dopamine activity^43^. This change was observed specifically within lNAc (the stimulation site) and not dorsal striatum (DS) (Fig. 2c), indicating that opto-conditioning was confined to the target region (lNAc). Selecting neurons inhibited by either odor did not reveal a clear enhancement of inhibition to the CS+; if anything, these neurons’ responses were also potentiated (Extended Data Fig. 4e), suggesting that the primary effect of opto-conditioning was to positively enhance odor-evoked spiking in a subset of neurons.

We next aimed to determine which striatal cell types contribute to the changes in striatal spiking generated by opto-conditioning. The two primary striatal cell types are D1- and D2-MSNs, expressing the D1 and D2 dopamine receptor, respectively. We used double transgenic mice to virally express a calcium indicator (GCaMP8m) in lNAc D1- or D2-MSNs along with ChrimsonR in VTA dopamine neurons (n=7 D1 mice, n=6 D2 mice; Fig. 2d, Extended Data Fig. 5a,b; Methods). These mice then underwent the same two-odor opto-conditioning paradigm as above. Consistent with asymmetric plasticity rules in D1- and D2-MSNs^44–46^, phasic dopamine stimulation selectively potentiated D1-MSN but not D2-MSN responses to preceding odor stimuli (Fig. 2e–i, Extended Data Fig. 5c–f). We conclude that the changes in lNAc spiking activity (Fig. 2b) primarily reflect changes in D1-MSNs.

A key feature of a value-like signal is that it reflects a prediction of upcoming reward but does not respond directly to reward itself. Consistent with this, D1-MSN odor responses gradually increased over training (Fig. 2h), but D1-MSNs did not respond directly to dopamine release (Fig. 2h, Extended Data Fig. 4c)^47^. We conclude that the D1-MSN response to the conditioned odor is “value-like” in that it reflects a gradually learned prediction of upcoming dopamine stimulation but does not respond to dopamine itself.

Together, these activity patterns are suggestive of a TD learning Dopamine release gradually potentiates the odor-evoked response in D1-MSNs (“value”); in turn, D1-MSNs stimulate lNAc-projecting VTA dopamine neurons via a temporal difference transformation, generating a positive dopamine response at onset and a negative omission dip at offset of their activity, closely recapitulating TD RPE (Fig. 2k).

Activation of lNAc D1-MSNs drives lNAc dopamine release via a TD transformation

To directly test this idea, we asked whether activation of D1-MSNs alone could generate the TD error-like patterns of dopamine activity observed in our opto-conditioning task, and more generally, whether D1-MSNs drive dopamine release according to a temporal difference transformation—the core computation in TD learning (red box, Fig. 1a). We prepared mice for simultaneous opto-stimulation of lNAc D1-MSNs with photometric recording of lNAc dopamine release and Neuropixels recordings of lNAc D1-MSNs (Fig. 3a) and designed 7 opto-stimulation patterns that generate both qualitative and quantitative predictions for a TD error-like response (Fig. 3b): (1) A 1 second pulse train at 20 Hz, (2) a 2 second pulse train at 20 Hz, (3) an upward ramp from 0 to 40 Hz, (4) a downward ramp from 40 to 0 Hz, and (5–7) 3 second pulse trains at 5, 10, and 20 Hz (all light pulses 1 ms at 40 mW/mm^2^).

The simulated TD error of these 7 patterns is shown in the right of Fig. 3b, with the temporal discount factor γ ranging from 1 (green; a pure derivative) to 0.1 (blue) in units of discount/sec. As has been previously noted^5,21,22^, TD error is similar (but not identical) to a time derivative. Indeed, in the schematics in Fig. 3b, the TD error shows several qualitative derivative-like features that are shared between TD in general (γ≤1) and the strict derivative (γ=1): (1) Transient positive/negative responses to the onset/offset of stimulation, which scale with the intensity of stimulation (patterns 5–7), (2) A prominent negative response to the offset of pattern 3, with a smaller positive response during the ramp-up phase, and (3) A prominent positive response to the onset of pattern 4, followed by a prolonged negative response as stimulation ramps down.

Remarkably, actual dopamine responses to D1-MSN stimulation displayed all of these derivative-like features (Fig. 3c, right). In contrast, D1-MSNs responded faithfully to the stimulation patterns themselves (Fig. 3c, left; Extended Data Fig. 6a). Notably, even though striatal MSNs are inhibitory neurons and directly project to VTA, the initial dopamine response was excitatory, especially for stimulation patterns with a rapid onset. This counterintuitive response pattern supports the existence of a derivative-like transformation between lNAc D1-MSNs and lNAc dopamine release, exactly as needed to implement TD learning (Fig. 1a).

While TD error is derivative-like, it deviates from a derivative when γ<1, most prominently in the steady-state response to sustained The steady-state response is zero for γ=1 (derivative) but negative for γ<1. Dopamine responses better matched TD, with negative steady state responses that scaled with stimulus intensity (Fig. 3c). Computationally, this key distinction gives the TD algorithm its preference for current over future rewards. These data raise the possibility that the temporal discount parameter γ could be set by a neural circuit mechanism between D1-MSNs and dopamine neurons—a possibility we explore in greater depth below.

Capturing the TD transformation as convolution with a biphasic linear filter

In theory, a TD transformation can be viewed as convolution with a biphasic filter, which applies positive weights to recent input, negative weights to preceding input, and sums these to generate the output (Fig. 3n)^25^. To describe our results quantitatively, we therefore used linear time-invariant (LTI) systems identification analysis, a standard approach in engineering^51^, to estimate the linear filter that best captures the transformation from an input (e.g. D1-MSN spiking) to an output (e.g. lNAc dopamine release) (Fig. 3d) (Methods).

Whereas the best-fit filter from stimulus to D1-MSN spiking resembled a delta function, scaling the input to produce the output (Extended Data Fig. 6d), the best-fit filter from D1-MSN spiking to dopamine release was biphasic, with a rapid positive component followed by a delayed negative component, thus implementing an approximate temporal difference (Fig. 3e; order 1 LTI model; Methods). The LTI models accurately predicted held-out data for both transformations (cross-validated R^2^: stimulation patterns → D1-MSNs = 0.95, D1-MSNs → dopamine release = 0.81). These results argue that our data can be adequately understood as the result of a linear system implementing a TD transformation between D1-MSNs and dopamine release.

To examine how well the dopamine response matches TD as opposed to more complex models, we fit LTI models of varying order, using a standard model formulation which assumes a particular functional form of the filter (Methods)^51^. Briefly, in this formulation the nth order LTI model is composed of n TD transformations, applied sequentially, followed by smoothing and scaling. Each TD transformation is given by convolution with the “TD filter,” hTD(t):

hTD(t)=δ˙(t)+lnγδ(t)

Where δ(t) and δ˙(t) are the Dirac delta function and its derivative, respectively. Convolving an input V(t) with this filter generates an output V˙(t)+lnγV(t), which is exactly the continuous form of the TD equation^21^. The TD LTI model is order It applies exactly 1 TD transformation, followed by smoothing and scaling. In contrast, the zeroth order model applies smoothing and scaling but lacks a TD transformation, whereas higher order models apply multiple TD transformations before smoothing and scaling, enabling higher order derivatives of the input to influence the output.

To capture the transformation from D1-MSNs to dopamine release, an order 1 LTI model was required, but higher order models did not improve performance (green lines, Fig. 3f). In contrast, a zeroth order LTI model was sufficient to capture the transformation from the stimulation patterns to D1-MSNs (red lines, Fig. 3f). This tightens the link between dopamine responses and the TD simpler models are insufficient, and more complex models are not needed.

Next, we used the LTI framework to demonstrate that dopamine responses are better fit by the TD equation than by a pure derivative. We refit the order 1 LTI model while fixing γ=1, thus implementing a pure derivative, as ln(1)=0. This model performed substantially worse on held-out data than the order 1 LTI model in which γ was free to vary (Fig. 3g).

The pattern of dopamine results was replicated in a larger cohort without D1-MSN recordings, receiving only patterns 1–4 (n=13 mice; Extended Data Fig. 6e,h,k), and no dopamine responses were observed in mice with tdTomato instead of ChrimsonR in D1-MSNs (n=6 mice; Extended Data Fig. 6f,i,l), ruling out optical artifacts or visual responses. Mice with GFP instead of GRABDA3m in VTA axons showed slow biphasic fluctuations in GFP fluorescence following opto-stimulation, potentially reflecting a hemodynamic response, but GFP signals were small compared to GRABDA3m (n=3 mice; Extended Data Fig. 6g,j,m).

TD-like response in optotagged lNAc-projecting VTA dopamine neurons

To gain insight into the TD computation at the level of single neurons in the VTA and better estimate the rapid dynamics of the dopamine response, we recorded spiking signals of antidromically optotagged lNAc-projecting VTA dopamine neurons while stimulating lNAc D1-MSNs with the same 7 stimulation patterns (Fig. 3h). The same TD-like patterns as in Fig. 3c were observed in dopamine spiking activity (Fig. 3i), well-captured by a biphasic filter with a sharper time course (cross-validated R^2^ = 0.69; Fig. 3i, right). Like photometry signals, single dopamine neurons required an order 1 LTI model, but not higher (Fig. 3j), and were better fit by TD than a pure derivative model (Fig. 3k). The latency of the onset of the excitatory response to pattern 4 was less than 50 ms, and the response peaked between 50 and 100 ms (Extended Data Fig. 6c). Importantly, these data localize the TD-like transformation downstream of D1-MSN spiking but upstream of VTA DA neuron spiking, and constrain the time course of activation, narrowing down potential mechanisms.

Extracting the temporal discount parameter γ from dopamine neuron responses to D1-MSN stimulation

The LTI modeling framework allows us to extract the temporal discount parameter γ directly from dopamine responses to D1 stimulation (Fig. 3l). We derived an exact mathematical relationship between γ and a parameter of the LTI model, namely γ=e-σ, where σ is the zero of the order 1 LTI model fit (Methods). Bootstrapped estimates of γ from optotagged dopamine neurons closely matched a previous estimate, measured in mice performing a cued delay task^48–50^ (γ^=0.53 discount/sec in our data, γ^=0.56 discount/sec in Masset et al.^48^; Fig. 3m; note, however, that in Sousa et al.^49^ γ^≈0.9). γ in the LTI model is linked to the ratio of the positive and negative lobes of the impulse response (P/N) via an exponential relationship (Fig. 3n,o). This points to the relative magnitude and timing of the net-excitatory and net-inhibitory components of the dopamine response to D1-MSN activation as a key determinant of γ, grounding a higher order cognitive parameter (the time horizon of reward preference) in a concrete neural circuit mechanism.

lNAc contributes a value-like signal to VTA dopamine

Because these experiments relied exclusively on artificial paradigms, we next asked whether lNAc contributes to dopamine TD RPE signaling in tasks with natural rewards. GCaMP6f was expressed transgenically in dopamine neurons by crossing DAT-Cre^+/−^ and Ai148^+/+^ mice. We then prepared mice for injection of muscimol (a GABAA receptor agonist) or saline into the lNAc with photometry recording in VTA (Fig. 4a). We chose to record somatic rather than axonal calcium signals because muscimol may directly inhibit dopamine axons, which express GABA receptors^52^.

Mice were water-deprived and trained in a classical conditioning task in which one odor (CS+) was associated with delivery of a natural reward (4μL water; delivered on 90% of CS+ trials) and another odor (CS−) was associated with nothing (Fig. 4b). In addition, unpredicted rewards were delivered without a preceding odor on ~5% of trials. After several days of training, we injected muscimol or saline into lNAc on alternating days while recording dopamine neuron calcium signals in the VTA. Muscimol dramatically increased calcium responses to predicted rewards, but not unpredicted rewards, and reduced responses to the CS+ (Fig. 4c–e, Extended Data Fig. 7). This pattern is consistent with a suppression of the reward expectation (value) associated with the CS+: the reward expectation of the CS+ decreases, so the cue response is suppressed, and at the same time the predicted (but not unpredicted) reward becomes more surprising, and so the reward response increases.

We tested for the specificity of this effect to odor value by also delivering novel odors in a subset of mice (4/6 mice). In contrast to the reduction in the dopamine response to the CS+, the response to novel odors did not change (Fig. 4e), arguing for a specific effect on odor value rather than novelty.

Because dopamine activity can reflect movements^53–55^, we analyzed face videos to check whether muscimol changed odor- and reward-triggered movements (Extended Data Fig. 8). Under muscimol, mice discriminated CS+ and CS− odors, consumed rewards, and detected novel odors (Extended Data Fig. 8b), suggesting that muscimol did not eliminate animals’ ability or motivation to perform the task. We did, however, detect changes in total mouth movement during the CS+ and after rewards and odor cues (Extended Data Fig. 8b), but muscimol similarly increased mouth movements to both predicted and unpredicted rewards, whereas VTA dopamine neuron calcium dramatically increased to predicted but not unpredicted rewards (Fig. 4e), indicating that movements cannot explain the main pattern of our results.

Together, these results demonstrate that lNAc contains a reward expectation (value)-like signal that contributes to TD RPE dopamine activity during classical conditioning with natural rewards, consistent with the TD error-like response of dopamine neurons to opto-stimulation of lNAc D1-MSNs.

The temporal difference computation is a widespread feature of striatal-dopamine circuitry

Dopamine neurons have heterogenous functions depending in part on their projection target^33,39,53,56–58^, and may not all signal TD errors or participate in value learning. To investigate whether the TD computation identified here generalizes to striatal dopamine systems beyond lNAc, we stimulated D1-MSNs in lNAc, dorsomedial striatum (DMS), and dorsolateral striatum (DLS) while recording dopamine release in each site (Fig. 5a,b, Extended Data Fig. 9). To our surprise, the estimated impulse response for each stimulation site was biphasic (Fig. 5c), indicating that the temporal difference transformation is a ubiquitous feature of the circuitry connecting striatal D1-MSNs to midbrain dopamine neurons. However, the parameters of the impulse response, including the estimated temporal discount factor γ, varied by striatal subregion, with the estimated γ being closer to 1 for DMS and DLS than lNAc (though not significantly so; Fig. 5d,e). This points to a unified computational account of phasic dopamine’s role in learning as a prediction error machine, with parameters and input/output spaces that vary by striatal subregion^59–61^.

Discussion

Collectively, our experiments reveal how key TD learning operations are accomplished by the circuitry connecting lNAc-projecting dopamine neurons and lNAc D1-MSNs: phasic lNAc dopamine release selectively potentiates lNAc D1-MSN responses to preceding sensory stimuli, and D1-MSNs drive VTA dopamine neurons spiking and lNAc dopamine release via a TD transformation, generating TD error-like activity (a “minimal TD learning loop”; Extended Data Fig. 10a). As a counterpoint to recent arguments over whether TD RL is a good description of dopamine’s role in learning at the behavioral level^15–18^, these results demonstrate that TD error calculations are hard-wired into the dopamine system at the microcircuit level.

While the response of dopamine neurons to D1-MSN stimulation is derivative-like, it is not a strict derivative, but instead more closely matches the temporal difference equation γVSt-VSt-1, with γ<1. This is revealed most clearly by the response to long pulse trains (patterns 5–7, Fig. 3c), to which the dopamine response plateaus at a negative value that becomes more negative as stimulation frequency increases. In the LTI framework, this corresponds to a biphasic impulse response with total negative area (N) greater than total positive area (P). The ratio of these areas (P/N) determines the temporal discount parameter γ of the TD algorithm (Fig. 3n,o).

A natural extension of these findings is that the same neural circuit mechanism underlies the cognitive phenomenon of temporal discounting—the tendency for animals and humans to prefer immediate over future rewards. The time horizon of this preference—set by γ—is a fundamental aspect of learning and decision making with implications for human health and economic choices^13,64–67^. For instance, discount rates for money correlate with individuals’ choices to engage in unhealthy behaviors, such as smoking^66^ and excessive alcohol consumption^68^. The biological basis of time discounting is not known. Our findings suggest a compelling γ is set at the level of individual dopamine neurons by the balance of the net-excitatory and net-inhibitory components of the circuit connecting D1-MSNs to VTA. This could explain why γ is a conserved property of single dopamine neurons across tasks^48^ and potentially explain the diversity of temporal discounting factors across striatal subregions^62,69^ and across individual animals in healthy and pathological conditions^66^. Notably, at the behavioral level temporal discounting is not exponential, as would be the case for a TD algorithm with a single discount factor, but hyperbolic, with a special status given to the current moment in time^12,50,70,71^ (but see^72^). A reservoir of dopamine neurons with private γ, each set individually by the balance the of positive and negative lobes of its impulse response to D1-MSN stimulation (Extended Data Fig. 10b), could combine to produce hyperbolic discounting on aggregate^48–50,73^. Further experiments linking the shape of the impulse response described here to neurophysiological and behavioral data in time discounting tasks would be needed to test this hypothesis (Extended Data Fig. 10c).

The TD transformation was a widespread feature of circuitry connecting D1-MSNs to dopamine release across both lNAc and dorsal striatum, revealing a conserved learning motif across striatal subsystems, despite their distinct roles in behavior^33,61^. This suggests that a single conserved TD-learning mechanism may have been repurposed by evolution to learn predictive maps across different input spaces^61,74^ and timescales^62,69^. Consistent with this, while the biphasic shape of the learned filter was conserved across the striatal areas we tested, the time resolution differed (time to minimum; Fig. 5e). The region-specific filters are likely shaped by a combination of both dopamine reuptake dynamics, which vary systematically across the striatum^63^, and circuit dynamics between striatum and dopamine neurons. Experiments involving projection-specific electrophysiological recording of the type developed here could distinguish between these alternatives.

Our experiments did not exhaustively cover striatal subregions, leaving several intriguing targets for future work, notably the medial nucleus accumbens, where dopamine release has value-like properties^39^, and the tail of the striatum, where dopamine may contribute to learning to predict threatening or salient outcomes^35,60,75^ or act as a value-neutral teaching signal^76^. If the hardwired TD computation extends to these regions, this would suggest dopaminergic learning follows an analogous algorithm despite notably different dopamine activity patterns. Outside striatum, the amygdala receives strong dopamine inputs from VTA, but appears to play a qualitatively distinct role in reinforcement learning compared to nucleus accumbens, with amygdala being more important for the initial acquisition of a behavioral response to motivationally salient cues (both rewarding and aversive)^77^, and nucleus accumbens subsequently refining the quantitative value estimate associated with these cues^78^. Distinct computations between amygdala/striatal neurons and dopamine neurons could help to explain this difference.

What are the possible biological mechanisms by which the temporal difference is computed? By recording lNAc D1-MSN optotagged neurons, antidromically optotagged lNAc-projecting VTA dopamine neurons, and lNAc dopamine photometry signals in the same paradigm, we show that the calculation occurs between the spiking of lNAc D1-MSNs and the spiking of lNAc-projecting VTA dopamine neurons, which argues against significant roles for local modulation of dopamine release^20,79^ or spiking adaptation in D1-MSNs in the TD error computation. We identify three categories of possible cell intrinsic mechanisms in dopamine neurons^80–82^, synaptic mechanisms in inputs to dopamine neurons^83,84^, and circuit mechanisms between D1-MSNs and dopamine neurons^27,83,85–87^. Cell intrinsic mechanisms could include depolarization block^82,88,89^ or D2 autoreceptor activation triggered by dopamine neuron spiking^90^. However, barring a spiking-independent mechanism, these mechanisms alone could not explain the pronounced suppression of action potential firing locked to the offset of pattern 3, which typically did not strongly activate spiking during the ramp-up phase (Fig. 3i). Next, synaptic mechanisms could include short-term synaptic depression (e.g. depletion of vesicle pools), which can create derivative-like responses in post-synaptic cells, wherein they respond proportionally to the percentage change of the input^82^. While short-term synaptic plasticity may play a role within a larger circuit computation^83^, a single depressing synapse cannot explain our finding that the dopamine response to lNAc D1-MSN stimulation plateaus at a level below baseline, which distinguishes TD error from a strict derivative. For these reasons, we hypothesize that circuit-level mechanisms play a major role in the temporal difference calculation. Importantly, although D1-MSNs are inhibitory neurons, activating them induced a burst of spikes in dopamine neurons (latency = 50–100 ms, Extended Data Fig. 6c), implying a disinhibitory mechanism. Disinhibitory circuit mechanisms could involve local inhibitory microcircuits within the VTA^83,85,91^ and neighboring substantia nigra pars reticulata (SNr), both of which receive innervation from lNAc D1-MSNs, or long-range multisynaptic basal ganglia circuits, such as from lNAc D1-MSNs to ventral pallidum (VP) to VTA/SNr^27,87^. Indeed, a recent study found that NAc D1-MSNs either (net) excite or (net) inhibit VTA dopamine neurons, depending on whether they project to the ventral midbrain (including VTA and SNr) or ventral pallidum^87^. If properly tuned in magnitude and timing, these oppositely signed responses could combine to form the temporal difference response reported here. The balance of these pathways would determine the temporal discount factor γ, thus grounding the TD learning algorithm’s sensitivity to time in a biological circuit mechanism (Extended Data Fig. 10b,c).

Altogether, our results show how TD computations are built into the circuitry of the basal ganglia, potentially explaining a wealth of neurophysiological observations, including ramping dopamine signals that possess a derivative-like property^5,19,20^. The same circuit mechanism may underlie temporal discounting at the level of single dopamine neurons^48,49^. While this shows that TD-like learning is a core property of the dopamine-striatum system when probed in isolation, natural behavior recruits additional circuits that may modulate or expand the functionality of the dopamine system. Therefore, the extent to which the LTI model proposed here can explain the contribution of dopamine and striatum to learning in more naturalistic contexts remains an important direction for future study.

Methods

Experimental procedures

Animals

A total of 74 mice (39 male, 35 female, aged 2–12 months at time of surgery) were used in the experiments. Mice were housed on a 12 hr / 12 hr dark/light cycle and experiments were performed in the dark phase. All procedures were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Harvard Institutional Animal Care and Use Committee.

Surgeries

All surgical procedures were performed in aseptic conditions. Mice were anesthetized using isoflurane (4% induction, 1–2% maintenance) and given an intraperitoneal injection of buprenorphine HCl (0.1 mg/kg) for analgesia. They were then head-fixed within a stereotactic surgical setup (Narishige). The hair was removed, the scalp was cleaned with 70% ethanol then 10% povidone-iodine, and a local anesthetic (2% lidocaine, 0.05 mL) was injected under the scalp. The scalp and periosteum were removed and the skull was dried and scored with a scalpel blade to ensure good adhesive binding. The skull was leveled by ensuring bregma and lambda were at the same DV position; no ML leveling was performed. Fiducial marks were made above the target sites (for virus injections, fiber implants, electrophysiology recording, or muscimol injection) using a fine-tipped pen. The following coordinates were used, with all distances in mm (AP and ML measured relative to bregma and DV measured from the surface of the brain): VTA, virus injection (Fig. 1–3, 5): −3.0 AP, 0.6 ML, 4.3→4.2→4.1 DVVTA, fiber implant (Fig. 4): −3.0 AP, 0.6 ML, 4.2 DVSNc, virus injection (Fig. 5): −3.0 AP, 1.6 ML, 4.2→4.1→4.0 DVlNAc virus injection (Fig. 1–3, 5): 1.3 AP, 1.8 ML, 4.1 DVlNAc fiber implant (Fig. 2, 3, 5): 1.3 AP, 1.8 ML, 4.0 DVlNAc angled fiber implant (Fig. 1–3): 1.3 AP, −1.1 ML, 4.75 DV, angled 35° in coronal plane away from midlineDMS, virus injection (Fig. 5): 0.5 AP, 0.5–0.7 ML, 2.7→2.6 DV, angled 15° in coronal plane towards midlineDMS, fiber implant (Fig. 5): 0.5 AP, 0.5–0.7 ML, 2.5 DV, angled 15° in coronal plane towards midlineDLS, virus injection (Fig. 5): −0.3 AP, 2.8 ML, 2.6 DV, angled 15° in sagittal plane towards posteriorDLS, mirrored fiber implant (Fig. 5): −0.3 AP, 2.6 AP, 2.8 DV, angled 15° in sagittal plane towards posterior, mirror facing towards lateral side

The exact combinations of viruses, fibers, and target sites depended on the experiment and are given in the corresponding sections. Craniotomies were made above virus injection and fiber implant target sites using a high-speed dental drill. Virus injections were performed using glass pipettes (Drummond, pulled using Narishige PC-10) with an oil hydraulic microinjector (Narishige MO-10). All virus injections were 300 nL volume with viral titer ranging from 5 × 10^12^ to 2 × 10^13^ gc/mL. The pipette was slowly withdrawn 10 minutes after completing the injection. Optic fiber cannulae were scored using the dental drill (to improve adhesive binding), slowly implanted using a stereotaxic cannula holder (SCH_1.25, Doric Lenses), and temporarily held in place using a small amount of superglue (n-butyl cyanoacrylate) and accelerant. After fiber implants, a custom titanium headplate was attached to the skull with dental cement (Metabond, Parkell) containing a small amount of charcoal powder to reduce light leakage. The cement also served to reinforce the optic fiber implants. For electrophysiology and muscimol experiments, a 3D-printed plastic well was glued to the headplate with epoxy or superglue, to hold saline during recordings or injections. For electrophysiology experiments, a ground pin coated with a lead-tin alloy (Mill-Max) was implanted above contralateral cortex. Apart from ground pin implants, all surgical procedures (virus injections, fiber implants, craniotomies) were done in the left side of the brain.

Following surgery, mice recovered on a heating pad, and were given ketoprofen (20 mg/kg, injected SC) and 0.5 mL of sterile saline (injected IP). Mice were then singly housed and post-operative analgesia (ketoprofen, 20 mg/kg) was given twice daily for 48 hours. In experiments using viruses, we waited at least 3 weeks for the virus to express before starting experiments.

Behavioral training

All behavioral experiments took place inside custom-built, partially sound-proofed, enclosed behavioral setups in which head-fixed mice ran freely on a cylindrical running wheel. White noise machines outside the setups were used to mask background noise. Except for opto-stim experiments using red masking lamps (Fig. 2j, Fig. 3, Fig. 5; see Red masking lamp below), experiments were run in the dark, with infrared lamps illuminating the mouse for face and body videos. Behavioral tasks were controlled using a Bpod state machine (1027, Sanworks) and accessories (1004: Port Interface, Optical Lickometer, Analog Output Module, Valve Module; Sanworks) with custom software written in MATLAB (R2021b) based on functionality in the Bpod MATLAB Software Repository (github.com/Sanworks/Bpod_Gen2).

Water deprivation

Mice were water-deprived starting at least 7 days after surgery. After this, mice were weighed daily and given water (1–2 mL) to maintain their bodyweight above 85% of baseline throughout the experiment.

Water delivery and lick detection

Water droplets were delivered to head-fixed mice through a spout attached to an optical lickometer (1020, Sanworks), which detected licks as breaks of an infrared beam. Water delivery was controlled by a Bpod port interface (1004, Sanworks) and valve (LHDA1231115H, The Lee Company). Water droplet volume was calibrated using a built-in Bpod protocol. In some experiments, to improve visibility of the face, the beam-break lickometer was removed and licks were detected directly from face videos in post-processing (e.g. Extended Data Fig. 8a, right).

Habituation to recording setups

All mice went through the same habituation protocol before starting 1 day of handling by experimenter, 1 day of freely exploring the head-fixation setup while receiving water from a syringe (administered by experimenter), 1 day of head-fixation while receiving water from a syringe (administered by experimenter), followed by several days of head-fixation with water droplets administered from the same valve/spout used in experiments. Thus, when mice began experiments, their weights had stabilized, and they were well habituated to handling by the experimenter, head-fixation, and receiving water droplets within the recording setup.

Odor preparation and delivery

Pure odors (1-hexanol, (S)-(−)-limonene, isoamyl acetate, p-cymene, ethyl butyrate, (+)-carvone, (+/−)-citronellal, eugenol, 1,4-cineole, L-fenchone, α-ionone; Sigma-Aldrich) were diluted to 10% in mineral oil (Sigma-Aldrich). Each day, 30μL of each odor solution used in the current experiment was placed on a syringe filter (Whatman Puradisc 6823–1327; 13 mm diameter, 2.7μm pore size) and the filters were attached to a custom odor manifold made of polyetheretherketone (PEEK). Filtered air (0.9 L/min) was constantly delivered to the mouse’s nose. A system of valves (LHDA1221111H, The Lee Company, mounted in a manifold, LFMX0510528B, The Lee Company), controlled by the Bpod state machine (1027, Sanworks) with an accessory valve module (1015, Sanworks), added 0.1 L/min to this air stream, for a total flow of 1 L/min. By default, this added air passed through a filter containing only mineral oil (“blank”). To deliver an odor, the valves switched the additional air flow from the blank filter to the desired odor filter. A vacuum line inside the recording setup provided air circulation to remove odors. In all experiments, odor identities assigned to CS+/CS− were counterbalanced across mice.

Opto-conditioning task

After standard habituation, mice started training on an odor-based opto-conditioning task (Fig. 1, Fig. 2). The two odors used were 1-hexanol and limonene. One odor was associated with opto-stimulation of dopamine axons in lNAc on 75% of trials (CS+) and the other was associated with nothing (CS−). Dopamine axon stimulation was calibrated to a natural reward (8μL water) for each mouse (Fig. 1e, Extended Data Fig. 1b, Extended Data Fig. 5b; see Fiber photometry with optogenetics for details). Odors were delivered for 1 second followed by a 0.5 second trace period (Fig. 1, Fig. 2) or 1 second trace period (Extended Data Fig. 1j–l) before opto-stimulation. Inter-trial intervals (ITIs) were drawn from a truncated exponential distribution with mean 4 seconds and maximum 12 seconds and added to a fixed 8.5 second interval for a total ITI of 8.5–20.5 seconds with mean 12.5 seconds. 80 trials of each type (CS+ and CS−) were delivered, pseudo-randomly interleaved, for a total of 60 CS+ stimulated trials, 20 CS+ omission trials, and 80 CS− trials. Immediately before and after the odor trials, we delivered blocks of water rewards, as a comparison point for the signals generated by odors and opto-stimulation. These blocks of water rewards consisted of 16–20 trials, split evenly between 2μL and 8μL droplets pseudo-randomly interleaved, with the same ITI distribution as the odor trials. At the very end of the session (after the second block of water rewards), 10 unpredicted opto-stimulation trials were delivered, with the same ITI distribution, to compare predicted and unpredicted opto-stimulation responses. Mice were trained in this task for up to 30 days, with photometry recording each day, and electrophysiology recordings on some days in late training (days 11+) (see Fiber photometry with optogenetics and Electrophysiology recording and spike sorting for details).

Classical conditioning task with water rewards

After standard habituation, mice started training on an odor-based classical conditioning task with water rewards (Fig. 4). The two odors used were isoamyl acetate and p-cymene. In the first phase, odors were paired with 100% reward delivery (4μL water) or nothing. Odors were delivered for 1 second followed by a 1 second trace period. Inter-trial intervals (ITIs) were drawn from a truncated exponential distribution with mean 4 seconds and maximum 12 seconds and added to a fixed 8.5 second interval for a total ITI of 8.5–20.5 seconds with mean 12.5 seconds. 80 trials of each type were delivered each day for 2–6 days. One mouse (MC190) started training with probabilistic reward delivery (80–90%) for 7 days before switching to 100% reward delivery for 3 days. After this shaping phase, mice were switched to the full task, in which the CS+ switched from 100% to 90% rewarded (in 1/10 trials, reward was omitted), and rare unpredicted rewards were delivered (1 in 21 trials). After training in the full task for 1–3 days, a craniotomy was performed over lNAc and muscimol or saline was injected into lNAc before the task each day for up to 8 days (see Muscimol injections). In some muscimol/saline sessions, we also delivered novel odors (5 or 10 trials total), either interleaved with other trial types, or at the end of the recording session. A different novel odor was used each day. The novel odors used 1-hexanol, (S)-(−)-limonene, ethyl butyrate, (+/−)-citronellal, eugenol, 1,4-cineole, L-fenchone, α-ionone.

Red masking lamp

In some experiments, including all experiments in which D1-MSNs were stimulated while recording dopamine signals (Fig. 3, Fig. 5) as well as the control experiment in Fig. 2j, the recording setups were illuminated with red light to mask visual input from red optogenetic stimulation light. Red light illumination came from two lamps with red light bulbs pointed directly at the mouse’s eyes. This eliminated visual dopamine responses to red light (Fig. 2j) so that optogenetic responses could be better isolated.

Face and body videos

In all experiments, face and body videos were recorded using infrared cameras (Blackfly BFLY-U3-03S2M, FLIR), running at least 30 frames per second, with custom software in Bonsai (https://bonsai-rx.org). Videos were synchronized to other data streams using TTL pulses output by BPod at the start of each trial, detected by cameras’ GPIO lines and saved to file for each video frame.

Fiber photometry with optogenetics

Photometry recordings were performed using a commercial photometry imaging system (Bundle-imaging Fiber Photometry System, Doric Lenses). A branching bundle patch cord (BBP(3)_400/440/900-0.37_1m_FCM-3xFC_LAF, Doric Lenses) was imaged by a CMOS camera (Sony IMX249LLJ) through an objective. The branching patch cord split into three branches, which could be used to image up to three sites simultaneously. A blue excitation LED (470 nm, 50μW power at tip of each 400μm diameter patch cord, 20 ms illumination at 20 Hz) was used to collect GRABDA3m (Fig. 1, 3, 5), GCaMP8m (Fig. 2), or GCaMP6f (Fig. 4) signals. A purple excitation LED (415 nm, 50μW power at tip of patch cord, 20 ms illumination at 20 Hz, interleaved with blue excitation pulses) was used to collect isosbestic signals, but these were only used to correct GCaMP signals as the isosbestic wavelength for GRABDA3m is not known. GRABDA3m recordings were not corrected by a control channel. Instead, control GFP mice were used to rule out movement artifacts, which were minimal in our head-fixed setup (Extended Data Fig. 1c). In experiments in which GRABDA3m and Neuropixels recordings were made at the same site (Fig. 2, Fig. 3), blue excitation light power was reduced to 20μW and purple excitation light was turned off to minimize light artifacts in the Neuropixels recording.

To perform photometry and optogenetics through the same optic fiber implant, red excitation light was independently coupled to each of the three patch cord branches using a fluorescence minicube (FMC3_(390–540)_(560–1000)_S, Doric Lenses). Each minicube contained a dichroic mirror which reflected blue/green photometry light and passed red optogenetic excitation light (cutoff wavelength ≈ 550 nm). Each minicube had three FC connectorized ports, which were connected (1) one branch of the bundle-branching patch cord, the other end of which connected to the Doric photometry setup, (2) a 635 nm red LED (CLED_635, Doric Lenses), and (3) a patch cord (MFP_400/440/110–0.37_2m_FC-ZH1.25(F)_LAF, Doric Lenses) connecting to the corresponding fiber implant on the animal. To further reduce optical crosstalk, the output of the red LED was filtered by a band-pass filter (Semrock FF01–630/92, transmitting 580–680 nm) installed in the corresponding minicube port. Each LED’s output was triggered by TTL pulses from the Bpod state machine, which controlled the behavioral task. With this setup, red optogenetic excitation light could be independently delivered to each of up to three optic fiber implants while simultaneously collecting photometry recordings. The lack of optical crosstalk was confirmed using recordings in tdTomato control animals (Extended Data Fig. 1b).

Electrophysiological recording

Spiking data were collected using Neuropixels 1.0 single shank probes^94^ or Neuropixels 2.0 four shank probes^95^. Neuropixels data were acquired via a National Instruments PXI chassis (NI PXIe-1061, with NI PXIe-8381 and NI PCIe-8381) with an IMEC PXIe card (IMEC PXIE_1000), recorded with SpikeGLX version 20190413-phase3B2 (https://billkarsh.github.io/SpikeGLX/). Craniotomies were performed the day before recording and covered with Kwik-Cast (World Precision Instruments), and mice were allowed to recover overnight. On the day of recording, a 3D-printed piece was placed over the head of the mouse, blocking from the mouse’s sight the experimenter manipulating the probe above the mouse’s head. Probes were coated in dye (DiD, DiI, or DiO; Thermo Fisher Scientific Vybrant Multicolor Cell-Labeling Kit) to visualize probe tracks in histology. Neuropixels probes were lowered using a Thorlabs micromanipulator (PT1-Z8) at 1–10μm/sec. Slower speeds were used for VTA insertions. The probe was inserted 100μm past the target depth, then retracted to the target depth, and allowed to settle for 10–30 minutes prior to starting the recording.

Optotagging dopamine neurons

Four shank Neuropixels 2.0 probes were used to record from antidromically optotagged VTA dopamine neurons projecting to lNAc (Fig. 1k–q, Fig. 3h–o). In the experiments of Fig. 1, ChrimsonR was expressed in dopamine neurons by injecting 300 nL of AAV8-hSyn-FLEX-ChrimsonR-tdTomato (UNC Vector Core) into the VTA of Dat-Cre^+/−^ mice. During the experiment dopamine axons in lNAc were illuminated with red light (635 nm LED, single pulses, 20 ms pulse duration, 40 mW/mm^2^) at the beginning and end of the recording session (40 trials each); these trials were used to identify lNAc-projecting dopamine neurons (see Definition of optotagged neurons). In the experiments of Fig. 3, ChR2 was expressed in dopamine neurons by injecting 300 nL of AAV5-Ef1a-fDIO-hChR2(H134R)-EYFP (UNC Vector Core) into the VTA of Dat-Flp^+/−^;Tac1Cre^+/−^ mice. Dopamine axons in lNAc were illuminated with blue light (473 nm laser, LaserGlow, single pulses, 20 ms pulse duration, 80 mW/mm^2^) at the beginning and end of the recording session (60 trials each).

Optotagging D1-MSNs

Single shank Neuropixels 1.0 probes were used to record from optotagged D1-MSNs in lNAc (Fig. 3a–g). ChrimsonR was expressed in D1-MSNs by injecting 300 nL of AAV8-hSyn-FLEX-ChrimsonR-tdTomato (UNC Vector Core) into the lNAc of Tac1-Cre^+/−^ mice. Optogenetic light pulses (635 nm LED, 40 mW/mm^2^) were 1 ms instead of 5 ms, to separate light artifacts from light-triggered spiking responses, which typically had latency ≈ 2–3 ms (Extended Data Fig. 6a). The first light pulse of the opto-stim patterns (Fig. 3a; excluding patterns 3 and 4) was used to identify optotagged D1-MSNs (see Definition of optotagged neurons).

Muscimol injections

Muscimol (Sigma-Aldrich M1523) was dissolved at 1ng/nL in sterile saline (0.9% NaCl) and aliquots were frozen at −20°C. Prior to experiments, an aliquot was thawed, further diluted to a final concentration of 0.5ng/nL in saline, and used for several days, stored at 4°C. One day before starting muscimol/saline injections, mice were anesthetized with isoflurane, and a craniotomy was made above lNAc (from AP: 1.3 mm, ML: 1.8 mm) and covered with KwikCast (World Precision Instruments). Mice were allowed to recover overnight. Injections of muscimol or saline (alternating days) then occurred the next day and every subsequent day for up to 8 days. For injections, mice were head-fixed and awake within the recording setup. The craniotomy was exposed, rinsed with saline, and visualized using a color camera (Chameleon3 CM3-U3-13S2C, FLIR). Muscimol or saline (50–100 nL) was injected into the lNAc (DV: 4.0 mm) using a glass pipette (Drummond, pulled using Narishige PC-10) and oil hydraulic microinjector (Narishige MO-10). After 10 minutes, the pipette was extracted from the brain. Mice then rested in their home cage for 1–2 hours before starting the behavioral task and photometry recording.

Multisite D1-MSN stimulation with GRABDA3m photometry

In this experiment (Fig. 5), a mirrored angled fiber was used to target DLS (MFC_400/430-0.48_3mm_MF1.25_MA45, Doric Lenses), facing towards the lateral side, to restrict red light illumination to DLS. Flat faced fibers were used to target lNAc and DMS. Mice had either 2 or 3 fibers implanted, targeting lNAc + DLS, or lNAc + DMS + DLS. An earlier cohort of mice had DLS implants which were not useable due to experimenter error (incorrect orientation of mirrored angled fiber); in these mice, only data from lNAc + DMS were analyzed.

Histology

Following experiments, mice were deeply anesthetized with an overdose of ketamine/medetomidine, exsanguinated with 0.9% phosphate buffered saline (PBS), and transcardially perfused with cold 4% paraformaldehyde (PFA) in PBS. The brain was extracted from the skull and stored in 4% PFA for 24–48 hours at 4°C, after which it was rinsed with PBS, stored in PBS, and cut into 100μm sections on a vibratome (VT1000S, Leica). In some experiments, sections were immunostained for tyrosine hydroxylase (AB152, Sigma Aldrich) to visualize dopamine axons or cell bodies. Sections were mounted onto glass slides with a mounting medium containing 4’,6-diamidino-2-phenylindole (DAPI; Vectashield, Vector Laboratories) and imaged with an automated fluorescence microscope (Axio Scan.Z1, Zeiss) or confocal microscope (LSM 700, Zeiss).

Data analysis

Unless otherwise noted, data were analyzed with custom code written in MATLAB (version 2024a).

Statistics

All statistical tests were performed in MATLAB (version 2024a). Unless otherwise noted, P values reflect two-tailed tests. We did not perform tests for normality or correct for multiple comparisons.

Data exclusion

In the muscimol experiment (Fig. 4, Extended Data Fig. 7, Extended Data Fig. 8), we excluded sessions in which injection (saline or muscimol) visibly failed (no movement of meniscus in glass pipette, possibly due to clogging). In total, 4/30 sessions were excluded (1 saline, 3 muscimol), leaving 26 sessions from 6 mice (16 saline, 10 muscimol). Including these sessions did not change the main pattern of results.

Photometry preprocessing

Raw photometry signals (F) were converted to deltaF/F (dF/F) through the formula dFF=F-F0F0, where F0 was defined as the 10^th^ percentile of F in a rolling 30 second window. dF/F traces were upsampled from 20 Hz to 1000 Hz through linear interpolation (MATLAB function interp1) and then smoothed with a Gaussian filter (SD = 50 ms). Smoothed dF/F traces were z-scored by subtracting the mean and dividing by the standard deviation of signals in the ITI periods, where ITI periods were defined by excluding time points from 0.5 seconds before to 3 seconds after the start of each trial. GCaMP signals (Fig. 2, Fig. 4) were corrected by regressing out the isosbestic signal, using linear regression (function robustfit in MATLAB) fit to ITI periods. GRABDA3m signals were not corrected by the isosbestic channel as the isosbestic wavelength for GRABDA3m is unknown, and using an incorrect isosbestic wavelength is more likely to distort the signal than correct it (motion artifacts were addressed for these animals using GFP control mice, e.g. Extended Data Fig. 1c and Extended Data Fig. 6g). Photometry data were synchronized to behavioral data using TTL sync pulses delivered to the photometry recording setup at the start of each behavioral trial.

Spike sorting

Neuropixels data were spike-sorted offline with Kilosort3^96^ using the ecephys pipeline (Allen Institute, adapted by Jennifer Colonell for SpikeGLX data; https://github.com/jenniferColonell/ecephys_spike_sorting), which includes CatGT to high-pass filter the raw data and detect and remove artifacts (https://github.com/billkarsh/CatGT). Kilosort clusters were then manually curated in Phy (https://github.com/cortex-lab/phy). After manual curation, quality control (QC) metrics and mean waveforms were computed for each unit using code from the ecephys pipeline. Neuropixels data were synchronized to behavioral data using TTL sync pulses delivered to the SMA input of the IMEC PXIe acquisition module at the beginning of each behavioral trial.

Peristimulus time histograms

Peristimulus time histograms (PSTHs) of neural activity (spike times or photometry) or behavioral data (motion energy) were computed using a MATLAB toolbox written by HyungGoo Kim (https://github.com/HRKimLab/libkm).

Definition of optotagged neurons

Optotagged dopamine neurons

Dopamine neurons were optotagged antidromically^97^ by stimulating dopamine axons in the striatum while recording spiking activity in the VTA with a 4 shank Neuropixels 2.0 probe (Fig. 1k,l). To be defined as optotagged, a unit had to meet the following Pass a stimulus-associated latency test (SALT)^98^ independently in two sets of 40 (Fig. 1k–q) or 60 (Fig. 3h–n) optotagging trials (single 20 ms pulses of light, red light for ChrimsonR in Fig. 1 and blue light for ChR2 in Fig. 3) delivered at the beginning and end of the recording session (SALT P<0.01 for both sets of trials). SALT compared spiking histograms from 1 to 20 ms (1 ms bins) following laser onset to a null distribution taken from ITI periods.Meet the following quality control criteria based on the QC metrics output by the ecephys pipeline, as recommended by the Allen Institute (https://allensdk.readthedocs.io/en/latest/_static/examples/nb/ecephys_quality_metrics.html): isi_viol < 0.5, amplitude_cutoff < 0.1, presence_ratio > 0.9. These metrics ensured neurons were well-isolated and present throughout the recording session.Have Pearson’s correlation at least 0.9 between mean waveform computed during light delivery (“light-triggered spikes”) and outside light delivery (“spontaneous spikes”). Waveforms were extracted using the software C_Waves written by Bill Karsh (https://github.com/billkarsh/C_Waves).Pass a collision test, a classic test of antidromic activation^97^. The first 1 ms bin following light onset in which firing rate differed significantly from a shuffle control was identified (z-score to shuffle > 5), which we may call tresp. Trials were identified in which spontaneous spikes occurred between tresp-7ms and tresp-2ms. To pass the collision test, these trials had to have 0 spikes in the 1 ms bin at tresp. If no spontaneous spikes occurred between tresp-7ms and tresp-2ms, the neuron was allowed to pass the collision test.Have mean firing rate > 0.2 Hz and < 10 Hz.

Optotagged D1-MSNs

D1-MSNs were optotagged by illuminating lNAc while recording spiking activity with a Neuropixels 1.0 probe (Fig. 3a–g). To be defined as optotagged, a neuron had to pass a stimulus-associated latency test (SALT)^98^ using the first pulse of each of the 7 stimulation patterns, excluding patterns 3 and 4 (SALT P<0.01). Light pulse duration was 1 ms in this experiment; SALT compared spiking histograms from 2 to 6 ms (1 ms bins) following laser onset to a null distribution taken from ITI periods.

Generating z-scored firing rate traces of single neuron spiking activity

Spike trains of individual neurons were transformed into smoothed z-scored firing rate traces by binning spikes (1 ms bins), smoothing the binned firing rate trace with a Gaussian kernel (SD = 50 ms), and z-scoring the resulting trace to ITI periods (Fig. 1o–q, Fig. 2b,c, Extended Data Fig. 2d–f, Fig. 3h–n, Extended Data Fig. 4). Here ITI was defined by excluding time points from 0.5 seconds before to 3 seconds after the start of each trial; all other time points were labeled “ITI”. Firing rate traces were z-scored by subtracting the mean and dividing by the standard deviation of ITI time points.

Definition of neurons significantly excited/inhibited by odors during opto-conditioning task

To define neurons as “significantly excited/inhibited” by an odor (CS+ or CS−) (Fig. 2b,c, Extended Data Fig. 4), trials were split into 1 second bins, spanning 1 second before odor onset to 5 seconds after odor onset. The number of spikes in each of the 5 1-second bins following odor onset was compared to the number of spikes in the 1-second bin preceding odor onset using a one-tailed ranksum test, separately for each odor (n = number of trials). Any neuron with P<0.001 in any of these 5 comparisons was deemed “significantly excited” (right tail; Fig. 2b,c, Extended Data Fig. 4a–d) or “significantly inhibited” (left tail; Extended Data Fig. 4e) by the odor.

Behavioral analysis

Video processing

Motion energy was extracted from within ROIs in each video. ROIs were drawn manually around the whisker pad, nose, mouth, and eye for face videos using custom MATLAB code. For body videos, the ROI was defined as the whole frame (full body motion). Within each ROI, motion energy was computed as the absolute value of the difference between subsequent frames.

Computational modeling

Linear time-invariant (LTI) systems identification analysis

Consider an input-output system ℱ, with input denoted by u(t) and output denoted by y(t), such that y(t)=ℱ(u(t)). This system is called linear time-invariant (LTI) if it has two Linearity:

(1) ℱ(a·u(t))=a·ℱ(u(t))

(2) ℱu1(t)+u2(t)=ℱu1(t)+ℱu2(t) Time

(3) ℱ(u(t-τ))=y(t-τ)

A system is LTI if and only if it is completely characterized by its impulse response, h(t), i.e. the response of the system to a Dirac delta function

(4) h(t)=ℱ(δ(t))

Where δ(t) is the Dirac delta function (not to be confused with the TD error, which we denote TD(t)). This means that (with zero initial condition), the output of an LTI system is equal to the convolution of the input with the impulse response (convolution denoted as *):

(5) y(t)=h(t)*u(t)

By taking the Laplace transform of both sides of equation 5, we can turn the convolution into a product (the convolution theorem of Laplace transforms):

(6) Y(s)=H(s)·U(s)

H(s), known as the transfer function of the system, is the Laplace transform of the system’s impulse response, and conversely, the impulse response h(t) is the inverse Laplace transform of the transfer

(7) H(s)=ℒ{h(t)}

(8) h(t)=ℒ-1{H(s)}

It is useful to consider the subset of all possible transfer functions H(s) which can be expressed as the ratio of two polynomials with real coefficients K, ai, and bj, with the numerator having order m and the denominator having order n:

(9) H(s)=Ksm+bm-1sm-1+⋯+b1s+b0sn+an-1sn-1+⋯+a1s+a0

This ratio-of-polynomials form can be factored as follows (fundamental theorem of algebra):

(10) =Ks-σ1s-σ2⋯s-σms-ρ1s-ρ2⋯s-ρn

Where these roots, σi and ρj, are complex numbers and are known as the “zeros” and “poles” of the system, respectively. As will be seen, the zeros and poles are very useful for interpreting and understanding the system dynamics. Many real-world LTI systems have H(s) with this form, but some do not (e.g. a strict time delay between input and output cannot be expressed this way).

When H(s) has this form, then the system is equivalent to a linear differential equation with constant coefficients, as

(11) Y(s)=H(s)·U(s)=Ksm+bm-1sm-1+⋯+b1s+b0sn+an-1sn-1+⋯+a1s+a0·U(s)

Rearranging

(12) snY(s)+an-1sn-1Y(s)+⋯+a0Y(s)=KsmU(s)+bm-1sm-1U(s)+⋯+b0U(s)

Taking the inverse Laplace transform of both

(13) y(n)(t)+an-1y(n-1)(t)+⋯+a0y(t)=Ku(m)(t)+bm-1u(m-1)(t)+⋯+b0u(t)

Where the last step follows from the derivative property of Laplace transforms ℒf′(t)=s·ℒ{f(t)}. Transforming the system into this form is useful for comparing it to the TD equation, which in the absence of rewards is a linear differential equation with m=1 (1 zero), n=0 (0 poles), K = 1, and b0=σ1=lnγ (see The TD equation is an LTI system in the absence of rewards).

Real biological systems are neither linear (e.g. neurons have threshold non-linearities) nor time-invariant (e.g. brain states like arousal fluctuate over seconds to minutes). We wished to understand the extent to which the transformation from striatal D1-MSN spiking to midbrain dopamine activity (spiking or release) can be described as an LTI system (i.e., how good of an approximation is this?), and to find the best fit h(t) describing this transformation (what is its shape? How many zeros and poles are needed to describe the system?). Finally, we wished to relate this empirical h(t) to the TD equation, which is itself an LTI system in the absence of rewards, transforming the value function V(t) into the TD error TD(t) by convolving with the TD impulse response, hTD(t) (see The TD equation is an LTI system in the absence of rewards).

The TD equation is an (improper) LTI system in the absence of rewards

The continuous form of the TD equation is^21^:

(14) TD(t)=r(t)+V˙(t)+lnγ·V(t)

In the absence of rewards (r(t)=0), this reduces

(15) TD(t)=V˙(t)+lnγ·V(t)

This is an improper (non-causal) LTI system with 0 poles and 1 zero, σ1=-ln(γ). The TD equation outputs zero

(16) V(t)=eσ1t=e(-lnγ)t=1γt

The TD equation with r(t)=0 is equivalent to convolving the value function V(t) with an impulse response function, which we call the “TD impulse response” or hTD(t):

(17) hTD(t)=δ′(t)+lnγ·δ(t)

(18) TD(t)=hTD(t)*V(t)

Where δ(t) is the Dirac delta function and δ′(t) is the derivative of the Dirac delta function.

Fitting LTI models to data

To fit LTI models to data (Fig. 3d–n, Fig. 5c–e, Extended Data Fig. 6e–m, Extended Data Fig. 9), we used MATLAB’s System Identification Toolbox, specifically the function tfest. LTI models were fit to the average PSTH for each stimulation pattern. To fit an order n model we used tfest with n zeros and n+1 poles (e.g. Fig. 3f). No regularization was used. Cross-validated R^2^ was computed by fitting the model on 6 of 7 stimulation patterns and predicting the response to the held-out pattern; these held-out predictions were concatenated across the 7 patterns, and cross-validated R^2^ was computed as 1-∑yi-yˆi2∑yi-y‾2, where yi is the true value of the average PSTH at time-point i, y^l is the prediction of the LTI model, and y‾ is the mean of all yi. To fit to a derivative model (order 1 model with σ1 fixed at 0), we first calculated the derivative of the input (MATLAB function diff), forcing σ1=0, and then fit an LTI model with 0 zeros and 2 poles. The cross-validated R^2^ of this model was compared against that of a full order 1 LTI model with 1 zero and 2 poles (Fig. 3g,k, Extended Data Fig. j-l).

Estimating γ from bootstrapped dopamine neuron population activity

We estimated γ from bootstrapped dopamine neuron population activity as follows. First, we sampled 55 optotagged dopamine neurons randomly with replacement from the full set of 55 neurons. Then we computed the average PSTH of this bootstrapped population and re-fit the LTI model (1 zero, 2 poles) predicting the average PSTH from the stimulation patterns. The estimated gamma, γˆ, was computed using the formula γˆ=e-σ1, where σ1 is the zero of the LTI fit (see above for derivation). This process was repeated 1000 times to generate the histogram in Fig. 3m. For each bootstrap, we also computed the ratio P/N (the ratio of the positive and negative areas of the fit impulse response). We plotted P/N against γˆ across all 1000 bootstraps in Fig. 3o. The black line in Fig. 3o shows an exponential curve of the form aebx+c, fit to the data using the MATLAB function fit.

Extended Data