Abstract
The neuromodulator dopamine plays a key role in motivation, reward-related learning, and normal motor function. The different affinity of striatal D1 and D2 dopamine receptor types has been argued to constrain the D1 and D2 signaling pathways to phasic and tonic dopamine signals, respectively. However, this view assumes that dopamine receptor kinetics are instantaneous so that the time courses of changes in dopamine concentration and changes in receptor occupation are basically identical. Here we developed a neurochemical model of dopamine receptor binding taking into account the different kinetics and abundance of D1 and D2 receptors in the striatum. Testing a large range of behaviorally-relevant dopamine signals, we found that the D1 and D2 dopamine receptor populations responded very similarly to tonic and phasic dopamine signals. Furthermore, because of slow unbinding rates, both receptor populations integrated dopamine signals over a timescale of minutes. Our model provides a description of how physiological dopamine signals translate into changes in dopamine receptor occupation in the striatum, and explains why dopamine ramps are an effective signal to occupy dopamine receptors. Overall, our model points to the importance of taking into account receptor kinetics for functional considerations of dopamine signaling.
SIGNIFICANCE STATEMENT Current models of basal ganglia function are often based on a distinction of two types of dopamine receptors, D1 and D2, with low and high affinity, respectively. Thereby, phasic dopamine signals are believed to mostly affect striatal neurons with D1 receptors, and tonic dopamine signals are believed to mostly affect striatal neurons with D2 receptors. This view does not take into account the rates for the binding and unbinding of dopamine to D1 and D2 receptors. By incorporating these kinetics into a computational model we show that D1 and D2 receptors both respond to phasic and tonic dopamine signals. This has implications for the processing of reward-related and motivational signals in the basal ganglia.
Introduction
The neuromodulator dopamine (DA) plays a key role in motivation, reward-related learning, and normal motor function. Many aspects of DA function are mediated by its effects on the excitability (Day et al., 2008) and strength of cortico-striatal inputs (Reynolds et al., 2001) in the context of motor control (Syed et al., 2016), action-selection (Redgrave et al., 2010), reinforcement learning (Schultz, 2007), and addiction (Everitt and Robbins, 2005). The striatal DA concentration ([DA]) can change over multiple timescales (Schultz, 2007). Fast increases in [DA] lasting for ≈ 1–3 s result from phasic bursts in DA neurons (Roitman et al., 2004), which signal reward-related information (Grace et al., 2007; Schultz, 2007). Slightly slower [DA] ramps occur as animals approach goal locations (Howe et al., 2013) or perform reinforcement learning tasks (Hamid et al., 2016). Finally, tonic firing of DA neurons may control the baseline [DA] and change on a timescale of minutes or longer (Grace et al., 2007). However, whether, e.g., learning and motivation are mediated by different timescales of DA cell firing (Niv et al., 2007) has recently been challenged (Berke, 2018; Mohebi et al., 2019). The issue of DA signaling time scales is important because D1 and D2 DA receptors may react to different timescales of the DA signal due to their different affinities for DA.
Based on the different DA affinities of the D1 receptor (D1R) and D2 receptor (D2R), it is often assumed that striatal medium spiny neurons (MSNs) respond differently to tonic and phasic DA signals, depending on which DA receptor type they predominantly express (Frank and O'Reilly, 2006; Grace et al., 2007; Schultz, 2007; Surmeier et al., 2007; Dreyer et al., 2010). According to this “affinity-based” model, the low-affinity D1Rs (high dissociation constant: KDD1 = 1.6μM; Richfield et al., 1989a) cannot detect tonic changes in [DA] because the fraction of occupied D1Rs is too small (≈1%) at a baseline [DA] of 20 nM and does not change much during tonic, low amplitude [DA] changes. However, D1Rs can detect phasic, high amplitude [DA] increases because they only saturate at a very high [DA]. By contrast, D2Rs have a high affinity (low dissociation constant: KDD2 = 25nM; Richfield et al., 1989a) leading to ≈40% of D2Rs being occupied at a baseline [DA] of 20 nM. Because of their high affinity, D2Rs can detect low amplitude, tonic [DA] increases/decreases. However, D2Rs saturate at relatively low [DA] > 2 · KDD2, and therefore cannot detect high amplitude, phasic [DA] increases. This suggests that D1 and D2 type MSNs respond differently to phasic and tonic [DA] changes because of the different affinities of D1Rs and D2Rs (Schultz, 2007). However, this model neglects other factors relevant for receptor occupation and is incompatible with recent findings that D2R-expressing MSNs can detect phasic [DA] changes (Marcott et al., 2014; Yapo et al., 2017).
The affinity-based model assumes that the reaction equilibrium is reached instantaneously, whereby the affinity can be used to approximate the fraction of receptors bound to DA. However, this assumption holds only if the receptor kinetics are fast compared with the timescale of the DA signal, which is typically not the case. For instance, D1Rs and D2Rs unbind from DA with a half-life time of t1/2 ≈ 80 s (Burt et al., 1976; Sano et al., 1979; Nishikori et al., 1980; Maeno, 1982), much longer than phasic signals of a few seconds (Robinson et al., 2001; Schultz, 2007; Hamid et al., 2016). Moreover, the fraction of bound receptors might be a misleading measure for the effect of DA signals, because the abundances of D1R and D2R in the striatum are quite different. Abundance here refers to the total number of D1 or D2 receptors available to bind to DA within a volume of extracellular space.
To investigate the role of receptor kinetics and abundances for DA signaling in the striatum, we developed a neurochemical model incorporating kinetics and abundances of D1Rs and D2Rs, and reevaluated current views on DA signaling in the striatum.
Materials and Methods
Code accessibility.
All models were implemented in Python. The models and all scripts used to generate the data and figures can be accessed here:
Kinetics model.
In the affinity-based model the receptor kinetics are instantaneous, so that the fraction of occupied D1 and D2 receptors (fD1 and fD2) can be calculated directly from the concentration of free DA in the extracellular space, [DA], and the dissociation constant KD (Copeland, 2004): However, the dissociation constant is an equilibrium constant, so it should only be used for calculating the receptor occupancy when the duration of the DA signal is longer than the time needed to reach the equilibrium. As this is typically not the case for phasic DA signals, because the half-life time of receptors is longer (Burt et al., 1976; Sano et al., 1979; Nishikori et al., 1980; Maeno, 1982) than the timeframe of phasic signaling (Roitman et al., 2004), we developed a model which incorporates slow kinetics.
When DA and one of its receptors are both present in a solution they constantly bind and unbind. During the binding process a receptor ligand complex (here called DA-D1 or DA-D2) is formed (Copeland, 2004). We refer to the receptor ligand complex as an occupied DA receptor. In the next section, refering to Eq. 2 up to Eq. 11, we provide the model equations for D1 receptors, but the same equations apply for D2 receptors (with different kinetic parameters). In a solution binding occurs when receptor and ligand meet due to diffusion, with high enough energy and a suitable orientation, described as follows: Accordingly, unbinding of the complex is denoted as follows: The kinetics of this binding and unbinding, treated here as first-order reactions, are governed by the rate constants kon and koff that are specific for a receptor ligand pair and temperature dependent. Because both processes are happening simultaneously we can write this as follows: The rate at which the receptor is occupied depends on [DA], the concentration of free receptor [D1] and the binding rate constant kon: The rate at which the receptor-ligand complex unbinds is given by the concentration of the complex [DA-D1] and the unbinding rate constant koff: The equilibrium is reached when the binding and unbinding rates are equal, so by combining Equations 5 and 6 we obtain the following: At the equilibrium the dissociation constant KD is defined as follows: When one-half of the receptors are occupied, i.e., [DA − D1] = [D1], Equation 8 simplifies to KD = [DA]. So at equilibrium, KD is the ligand concentration at which one-half of the receptors are occupied.
Importantly, for fast changes in [DA] (i.e., over seconds) it takes some time until the changed binding (Eq. 5) and unbinding rates (Eq. 6) are balanced, so the new equilibrium will not be reached instantly. The timescale in which equilibrium is reached can be estimated from the half-life time of the bound receptor. The half-life time assumes an exponential decay process as described in Equation 6 and is the time required so that one-half of the currently bound receptors unbind. Although the t1/2 estimates have been published, it is also possible to calculate it from experimental estimates of koff by using t1/2 = ln(2)/koff (Burt et al., 1976; Sano et al., 1979; Nishikori et al., 1980; Maeno, 1982). Signal duration should be of the same order of magnitude or longer than the half-life time for the affinity-based model with instant kinetics to be applicable.
We calculated the time course of occupied receptor after an abrupt change in [DA] by integrating the rate equation, given by the sum of Equations 5 and 6: To integrate Equation 9 we substitute the following: where [D1tot] is the total amount of D1 receptor (bound and unbound to DA) on the cell membranes available for binding to extracellular DA.
To model the effect of phasic changes in [DA] we choose the initial receptor occupancy [DA − D1](t = 0) = [DA − D1]0 and the receptor occupancy for the new equilibrium at time infinity [DA − D1](t = ∞) = [DA − D1]∞ as the boundary conditions. With these boundary conditions we get an analytic expression for the time evolution of the receptor occupancy under the assumption that binding to the receptor does not significantly change the free [DA]: For arbitrary DA time courses we solve Equation 9 for each receptor type numerically using a fourth-order Runge–Kutta solver with a 1 ms time resolution.
We did not take into account the change in [DA] caused by the binding and unbinding to the receptors because the rates at which DA is removed from the system by binding to the receptors is much slower than the rate of DA being removed from the system by uptake through DA transporters. The rate at which DA binds to the receptors is as follows: Equation 12 relates the removal of DA to the binding of DA to its receptors. To estimate the binding of DA to its receptors we use the parameters of the DA step example (Fig. 1). In this example there is an instantaneous DA increase from the baseline value [DA] = 20nM to [DA] = 1 μM. At the time of the step up, the D1 and D2 occupancy is given as [DA − D1] ≈ 20.0 nM and [DA − D2] ≈ 40 nM (the equilibrium values for baseline DA). With that the free receptor concentration is [D1] = [D1]tot − [D1 − DA] ≈ 1600.0nM and [D2] = [D2]tot − [D2 − DA] ≈ 40.0nM. The receptor parameters konD1 = 5.2 · 10−6nM−1s−1, konD2 = 3.3 · 10−4nM−1s−1, and the receptor abundances [D1]tot and [D2]tot are derived in the Receptor parameters section. For the parameters from this example, the rate of DA removal through binding to the receptors is given by the following: However, the DA removal rate by Michaelis-Menten uptake through the DA transporters at this concentration would be: Vmax = is the maximal uptake rate, and Km = 0.21 μM the Michaelis–Menten constant describing the [DA] concentration at which uptake is at half of the maximum rate (Bergstrom and Garris, 2003). Because , the DA dynamics are dominated by the uptake process and not by binding to the receptors. Therefore, we neglected the receptor-ligand binding for the DA dynamics in our model. However, for faster DA receptors this effect would become more important.
Receptor parameters.
The DA receptor abundances, i.e., the total concentration of the D1 and D2 receptors on the membrane ([D1]tot and [D2]tot) that can bind to DA in the extracellular space of the striatum are important model parameters. Our estimate of [D1]tot and [D2]tot is based on radioligand binding studies in the rat rostral striatum (Richfield et al., 1987, 1989a). We use the following equation in which X is a placeholder for the respective receptor type, to calculate these concentrations. The experimental measurements provide us with the number of receptors per unit of protein weight [D1]m and [D2]m. To transform these measurements into molar concentrations for our simulations, we multiply by the protein content of the wet weight of the rat caudate nucleus ϵ, which is ∼12% (Banay-Schwartz et al., 1992). This leaves us with the amount of protein per gram of wet weight of the rat brain. Next we divide by the average density of a rat brain, which is ρb = 1.05 g/ml (DiResta et al., 1990) to find the amount of receptors per unit of volume of the rat striatum. Finally, we divide by the volume fraction α, the fraction of the brain volume that is taken up by the extracellular space in the rat brain, to obtain the receptor concentration of the receptor in the extracellular medium. The procedure ends here for the D1 receptors because there is no evidence that D1 receptors are internalized in the baseline state (Nishikori et al., 1980; Prou et al., 2001). However, a large fraction of the D2 receptors is retained in the endoplasmic reticulum of the neuron (Nishikori et al., 1980; Prou et al., 2001), reducing the amount of receptors that contribute to the concentration of receptors in the extracellular medium by fm, the fraction of receptors protruding into the extracellular medium. Thus, we define the receptor abundances [DX]tot, the total number of receptors per unit of volume in the extracellular medium for both receptor types. It is useful to give this quantity as a concentration (nM), so that it can be easily used for calculating binding rates and equilibria as described in the beginning of this section.
Nishikori et al. (1980) also give estimates corresponding to [D1]m and [D2]m. Their measurements give [D1]m = 11.9 pmol/mg (protein) and [D2]m = 0.16 pmol/mg (protein). Because these measurements are already for the cellular membrane (i.e., they have already been implicitly multiplied by fmembrane), we used Equation 17 to obtain [D1]tot ≈ 6400 nM and [D2]tot ≈ 100 nM. Note that estimates of [D1]m may differ by a factor three among in the nucleus accumbens of rats, cats, and monkeys (Richfield et al., 1987). Here, we used the values corresponding to rats (Richfield et al., 1987, 1989b) instead of canine (Nishikori et al., 1980). Despite the differences across species, the receptor abundances derived from Nishikori et al. (1980) indicate that, at baseline [DA], the [D1–DA] would still be of the same order of magnitude as [D2–DA]. However, with these values [D1–DA] would be ∼twice [D2–DA] (as opposed to one-half in our model; Fig. 1).
In addition to the receptor concentration, the kinetic constants of the receptors are key parameters in our slow kinetics model. In an equilibrium measurement in the canine caudate nucleus the dissociation constant of low affinity DA binding sites, corresponding to D1 receptors (Maeno, 1982), has been measured as KD = 1.6μM (Sano et al., 1979). However, when calculating KD (using Eq. 8) from the measured kinetic constants (Sano et al., 1979) the value is KDD1 = 2.6 μM. To be more easily comparable to other simulation works (Dreyer et al., 2010) and direct measurements (Sano et al., 1979; Richfield et al., 1989a) we choose KDD1 = 1.6 μM in our simulations. For this purpose we modified both the konD1 = 0.00025min−1nM−1 and koffD1 = 0.64min−1 rate measured (Sano et al., 1979) by ≈25%, making konD1 = 0.0003125min−1nM−1 slightly faster and koffD1 = 0.5 min−1 slightly slower, so that the resulting KDD1 = 1.6 μM. The kinetic constants have been measured at 30°C and are temperature dependent. In biological reactions a temperature change of 10°C is usually associated with a change in reaction rate around a factor of 2–3 (Reyes et al., 2008). However, the conclusions of this paper do not change for an increase in reaction rates by a factor of 2–3 (see Fig. 9). It should also be noted that the measurements of the commonly referenced KD (Richfield et al., 1989a) have been performed at room temperature.
The kinetic constants for the D2 receptors were obtained from measurements at 37°C of high affinity DA binding sites (Burt et al., 1976), which correspond to the D2 receptor (Maeno, 1982). The values are konD2 = 0.02min−1nM−1 and koffD2 = 0.5min−1, which yields KDD2 = 25nM, in line with the values measured by Richfield et al. (1989a). As the measured off-rate of the D1 and D2 receptors koffD1 = 0.64min−1 and koffD2 = 0.5min−1 is quite similar and we modify the measured values slightly in our model (see previous paragraph), the difference in KDD2 = 25nM and KDD1 = 1.6 μM is largely because of differences in the on-rate of the receptors. This is important because the absolute rate of receptor occupancy depends linearly not only on the on-rate, but also on the receptor concentration (Eq. 5), which means that a slower on-rate could be compensated for by a higher number of receptors.
The parameters used in the simulations are summarized in Table 1.
Dopamine signals.
In our model we assumed a baseline [DA] of [DA]tonic (Suaud-Chagny et al., 1992; Justice, 1993; Venton et al., 2003; Borland et al., 2005; Dreyer et al., 2010; Dreyer, 2014; Atcherley et al., 2015). We modeled changes in [DA] to mimic DA signals observed in experimental studies. We use three types of single pulse DA signals: (long-)bursts, burst-pauses, and ramps.
The (long-)burst signal mimics the effect of a phasic burst in the activity of DA neurons in the SNc, e.g., in response to reward-predicting cues (Pan et al., 2005). The model burst signal consists of a rapid linear [DA] increase (with an amplitude ΔDA and rise time trise) and a subsequent return to baseline. The return to baseline is governed by Michaelis Menten kinetics with appropriate parameters for the dorsal striatum Vmax = 4.0 μMs−1 and Km = 0.21 μM (Bergstrom and Garris, 2003) and the nucleus accumbens Vmax = 1.5 μMs−1 (Dreyer and Hounsgaard, 2013). In our model the removal of DA is assumed to happen without further DA influx into the system (baseline firing resumes when [DA] has returned to its baseline value). Unless stated otherwise, the long-burst signals are used with a Δ[DA] = 200 nM and a rise time of trise = 0.2 s at Vmax = 1.5 μMs−1, similar to biologically realistic transient signals (Robinson et al., 2001; Cheer et al., 2007; Day et al., 2007).
The burst-pause signal has two components, an initial short, small amplitude burst (Δ[DA] = 100nM, trise = 0.1s), with the corresponding [DA] returning then to baseline (as for the long burst above). However, there is a second component in the DA signal, in which [DA] falls below baseline, simulating the effect of a pause in DA neuron firing. The length of this firing pause is characterized by the parameter tpause. We simulated this type of burst-pause [DA] signal to investigate how, e.g., a two-component response of the DA reward prediction error (Schultz, 2016) would affect the DA receptor occupation. In this case the model input [DA] time course is based on DA cell firing patterns consisting of a brief burst followed by a pause in activity (Pan et al., 2008; Schultz, 2016).
The ramp DA signal is characterized by the same parameters as the burst pattern, but with a longer trise, and a smaller Δ[DA] (parameter settings provided in each simulation).
For the simulations comparing the area under the curve of the input DA signal with the resulting receptor occupancy (see Fig. 5) we used the burst, burst-pause, and ramp signals described with a range of parameter settings. For the burst DA signal we used amplitudes Δ[DA]max of 50, 100, 150, 200, 250, 300, 350, 400, 500, 600, 700, 800, 900, and 1000 nM. For the ramping DA signals we used rise times trise of 0.2, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, and 7.0 s, For the burst-pause DA signal we used different values for Vmax of 1.0, 1.5, 2.0, 2.5, 3.0, 2.5, and 4.0 μM · s−1.
Behavioral task simulation.
To determine whether DA receptor occupancy can integrate reward signals over minutes, we simulated experiments consisting of a sequence of 50 trials. In each sequence the reward probability was fixed. The trials contained either a (long-)burst DA signal (mimicking a reward) or a burst-pause DA signal (mimicking no reward) at the beginning of the trial according to the reward probability of the sequence. The intertrial interval was 15 ± 5 s (see Fig. 8). We choose this highly simplistic scenario to mimic DA signals in a behavioral task in which the animal receives unpredictable rewards with a given reward probability. Because of the unpredictable nature of the reward, we assumed that here the DA pulse amplitude is not affected by the reward probability, because, e.g., DA cell recordings during un-signaled reward presentations on a similar time scale have been used to obtain strong DA cell responses (Fiorillo et al., 2008).
However, here the specifics of the task are not relevant as our model addresses the integration of the DA receptor occupancy over time. Although we chose to use the burst-pause type signal as shown in Figure 2a as a non-rewarding event, the difference to a non-signal are minimal after the end of the pause (Figs. 3, 4). Each sequence started from a baseline receptor occupancy, assuming a break between sequences long enough for the receptors to return to baseline occupancy (∼5 min). For the simulations shown in Figure 4 all trials started exactly 15 s apart.
Although for the simulations shown in Figure 4 the sequence of DA signals was fixed, we also simulated a behavioral task with stochastic rewards (see Fig. 8). There we simulated reward probabilities from 0 to 100% in 10% steps. For each reward probability we ran 500 sequences, and calculated the mean receptor occupancy over time (single realizations shown in Fig. 8a,c). To investigate whether the receptor occupancy distinguished between different reward probabilities we applied a simple classifier to the receptor occupancy time course.
The classifier was used to compare two different reward probabilities at a time. At each time point during the simulated experiment it was applied to a pair of receptor occupancies, e.g., one belonging to a 70% and one to a 30% reward probability sequence. For each sequence the classifier assigned the current receptor occupancy to the higher or lower reward probability depending on which reward probabilities' mean receptor occupancy (>500 sequences) was closer to the current receptor occupancy.
Because we knew the underlying reward probability of each sequence we were able to calculate the true and false-positive rates for each time point in our set of 500 sequences for both the D1R and D2R (see Fig. 8e). For each individual comparison of two reward probability sequences, at each time point the classifier could make a correct or incorrect classification (denoted by a “1” or “0”, respectively). The true and false-positive rates were then obtained by averaging these sequences of zeros and ones at each time point over 500 different realizations of the two reward probability sequences. For example, Figure 8e shows that at a time of 400 s in ∼94% of the 1000 (500 for each probability) studied sequences the underlying reward probability was correctly identified by the classifier (i.e., 30% as 30%, and 70% as 70%) based on the D1 receptor occupation. The classification was done by calculating the distance of the instantaneous D1 receptor occupancy (i.e., in this case at 400 s) of a given individual sequence to the mean (over 500 sequences) receptor occupation of the 30 and 70% cases and then choosing the closer one. From the time resolved true and false-positive rates we calculated the time averaged true and false-positive rates, for all pairs of probabilities (see Fig. 8b,d) and the time averaged accuracy (8f) using all time points between 200 and 800 s within a sequence to avoid the effect of the initial “swing-in” and post-sequence DA levels returning to baseline.
Results
Before investigating the role of the receptor kinetics in response to different DA signals, we started by establishing the receptor binding at baseline [DA], taking into account the different abundances of D1 and D2 receptors in the striatum. For a stable baseline [DA] the resulting receptor occupation can be calculated using the receptor affinities (see Materials and Methods; Eq. 1). We report the resulting receptor occupation as the concentration of D1Rs and D2Rs bound to DA (denoted as [D1–DA] and [D2–DA], respectively). Expressing receptor occupation in terms of concentration (typically in nM) follows from our estimates based on experimental measurements (Richfield et al., 1987, 1989a), and is convenient for the calculation of binding rates and equilibria (see Materials and Methods).
First, we investigated receptor binding for a range of affinities (Fig. 1), reflecting the range of measured values in different experimental studies (Neve and Neve, 1997). Because of the low affinity of D1Rs, at low baseline [DA] only a small fraction of D1 receptors may be occupied. However, there are overall more D1Rs than D2Rs (Richfield et al., 1989a), and ≈80% of D2Rs are retained in the endoplasmic reticulum (Prou et al., 2001). Therefore, the concentration of D1Rs in the membrane available to extracellular DA is a lot higher than the concentration of D2Rs (e.g., 20 times more in the nucleus accumbens; Nishikori et al., 1980; see Materials and Methods). Thus, in our simulation, the actual concentration of bound D1Rs ([D1 − DA] ≈ 20 nM) was, at DA baseline, much closer to the concentration of bound D2Rs ([D2 − DA] ≈ 35 nM) than suggested by the different D1 and D2 affinities alone. We further confirmed that this was not due to a specific choice of the dissociation constants in the model, as [D1–DA] and [D2–DA] remained similar over the range of experimentally measured D1R and D2R affinities (Neve and Neve, 1997); Fig. 1a. This suggests that [D1–DA] is at most twice as high as [D2–DA] instead of 40 times higher as suggested by the difference in fraction of bound receptors. Therefore, [D1–DA] and [D2–DA] might be better indicators for the signal transmitted to MSNs, as the fraction of bound receptors neglects the different receptor type abundances.
Next, we investigated the effect of slow [DA] changes (Grace, 1995; Schultz, 1998; Floresco et al., 2003) by exposing our model to changes in the [DA] baseline. For signaling timescales that are long with respect to the half-life time of the receptors (tslow ≫ t1/2 ≈ 80 s), we used the dissociation constant to calculate the steady-state receptor occupancy. We found that for a range of [DA] baselines (mimicking slow changes in [DA]), there was less than a twofold difference between [D1–DA] and [D2–DA] (Fig. 1b), because of the different abundances of D1 and D2 receptors. This is in contrast to affinity-based models, which suggest that D2Rs are better suited to encode slow or tonic changes in [DA]. Interestingly the change of [D1–DA] was almost linear in [DA], while the change of [D2–DA] showed nonlinear effects due to the change in available free D2R. Thus, based on these results, it could even be argued that D1Rs are better at detecting tonic signals at high [DA] levels, because they do not saturate as easily.
Although for baseline and slow changes in [DA] the receptor occupation can be determined based on the receptor affinity, fast changes in [DA] also require a description of the underlying receptor kinetics. To investigate the effect of typical DA signals on receptor occupation, we developed a kinetics model incorporating binding and unbinding rates that determine the overall receptor affinity (see Materials and Methods; Eqs. 8, 9). The available experimental measurements indicate that the different D1R and D2R affinities are largely because of different binding rates, while their unbinding rates are similar (Burt et al., 1976; Sano et al., 1979; Maeno, 1982; Richfield et al., 1989a). We incorporated these measurements into our kinetics model and investigated the model's response to a variety of fast DA signals.
We started by measuring the model response to a [DA] step change from 20 nm to 1 μM. This is quite a large change compared with phasic DA signals in vivo (Robinson et al., 2001; Cheer et al., 2007; Hamid et al., 2016), which we choose to illustrate that our results are not just due to a small amplitude DA signal. We found that binding to both receptor subtypes increased very slowly. Even for the high affinity D2Rs it took >5 s to reach their new equilibrium (Fig. 1c). Thus, unlike the affinity-based model, our model suggests that the D2Rs will not saturate for single reward events, which last overall for up to ≈3 s. Note that the non-saturation is independent of the abundance of the receptors and is only determined by the kinetics of the receptors (see Materials and Methods). Due to their slow unbinding, D1Rs and D2Rs also took a long time to return to baseline receptor occupancy after a step down from [DA] = 1 μM to [DA] = 20 nM (Fig. 1d). Thus, we conclude that with slow kinetics of receptor binding both D1Rs and D2Rs can detect single phasic DA signals and that both remain occupied long after a high [DA] has returned to baseline.
DA receptor binding kinetics for different types of DA signals
Next, we investigated [D1–DA] and [D2–DA] for three different types of DA signals (Fig. 2). The first signal was a phasic DA increase (Fig. 2a, long burst), mimicking responses to rewards and reward-predicting stimuli (Robinson et al., 2001; Cheer et al., 2007). The second signal was a brief phasic DA increase, followed by a decrease (Fig. 2a, burst-pause), mimicking responses to conditioned stimuli during extinction (Pan et al., 2008) or to other salient stimuli (Schultz, 2016). The third signal was a prolonged DA ramp, mimicking a value signal when approaching a goal (Howe et al., 2013; Hamid et al., 2016; Fig. 2b). In the affinity-based model with instant kinetics the D1Rs mirrored the [DA] time course for all three types of signals, because even at [DA] = 200 nM D1Rs are far from saturation. By contrast, D2Rs showed saturation effects as soon as [DA] > 2 · KDD2, leading to differing D1 and D2 time courses (Fig. 2, gray traces). Importantly, in our model with slow kinetics, the time courses of [D1–DA] and [D2–DA] were nearly identical (Fig. 2, bottom row), supporting that both receptor types are equally affected by phasic DA signals. This was the case for all the three signals: burst, burst-pause and ramping DA signals. The only difference between the [D1–DA] and [D2–DA] time courses were the absolute amplitudes. For example, [D2–DA] started from a baseline ∼twice as high as [D1–DA], but then also responded to the long burst DA signal with a change ∼twice as high. The similarity of [D1–DA] and [D2–DA] responses to both slow (Fig. 1b) and fast (Fig. 2) [DA] changes indicates that the different DA receptor types respond similarly independent of the timescale of [DA] changes. It could even be argued that D2Rs are better at detecting phasic DA signals, because they respond with a larger absolute change in occupied receptors.
To understand why the D1Rs and D2Rs respond in a similar fashion, we considered the relevant model parameters in more detail. The binding rate constants of D1Rs and D2Rs differ by a factor of ≈60 (konD1 = 0.0003125nM−1min−1 and konD2 = 0.02 nM−1min−1; Burt et al., 1976; Sano et al., 1979; Maeno, 1982; see Materials and Methods), suggesting faster D2Rs. However, experimental data indicates that there are ≈40-fold more unoccupied D1 receptors ([D1] ≈ 1600nM) than unoccupied D2 receptors ([D2] ≈ 40nM) on MSN membranes in the extracellular space of the rat striatum. This difference is due to a combination of simply higher abundance of D1 receptors (Richfield et al., 1987, 1989b) and D2 receptors being retained in the endoplasmic reticulum (Prou et al., 2001). Therefore, the absolute binding rate = kon · [DA] · [DX] differs only by a factor of ≈1.5 between the D1Rs and D2Rs. That is, the difference in the kinetics of D1Rs and D2Rs is compensated by the different receptor numbers, resulting in nearly indistinguishable aggregate kinetics (Fig. 2). This is consistent with recent experimental findings that D2R-expressing MSNs can detect phasic [DA] signals (Marcott et al., 2014; Yapo et al., 2017).
The dynamics introduced by the slow kinetics in our model also affected the time course of DA signaling. With instant kinetics the maximum receptor occupancy was reached at the peak [DA] (Fig. 2). By contrast, for slow kinetics the maximum receptor occupancy was reached when [DA] returned to its baseline (Fig. 2a) because as long as [DA] was higher than the equilibrium value of [D1–DA] and [D2–DA], more receptors continued to become occupied. Therefore, for all DA signals, the maximum receptor occupancy was reached toward the end of the pulse.
Another striking effect of incorporating receptor kinetics was that a phasic increase in [DA] kept the receptors occupied for a long time (Fig. 2a, green traces). However, when a phasic increase was followed by a decrease, [D1–DA] and [D2–DA] returned to baseline much faster (Fig. 2a, purple traces). This indicates that burst-pause firing patterns can be distinguished from pure burst firing patterns on the level of the MSN DA receptor occupancy. Furthermore, this supports the view that the fast component of the DA firing patterns (Schultz, 2016) is a salience response, and points to the intriguing possibility that the pause following the burst can, at least partly, revoke the receptor-ligand binding induced by the burst. In fact, for each given burst amplitude, a sufficiently long pause duration could cancel the receptor activation (Fig. 3), with larger [DA] amplitudes requiring longer pauses to cancel the activation. Thereby, the burst-pause firing pattern of DA neurons could effectively signal a reward “false-alarm”.
The long time it took [D1–DA] and [D2–DA] to return to baseline after phase DA signals (Fig. 2a) indicates that the receptor occupation integrates DA signals over time. To examine this property, we simulated a sequence of DA signals on a timescale relevant for behavioral experiments (Fig. 4). Each sequence consisted of 50 events and each event was separated by 15 s. Three different types of sequences were tested: 50 phasic DA bursts, 40 phasic DA bursts followed by 10 burst-pause signals, and 40 phasic DA bursts followed by 10 non-events. We found that both [D1–DA] and [D2–DA] accumulated over the sequence of DA signals. The sawtooth shape of the curves was due to the initial unbinding of the receptors after each burst event, which was then interrupted by the next DA signal 15 s later. At higher levels of receptor activation, the amount of additional activated receptor per DA pulse was reduced because there are less free receptors available, and the amount of receptors unbinding during the pulse duration was increased because more receptors were occupied. The accumulation occurred as long as the time interval between the DA signals was shorter than ≈ 2 · t1/2. Together, the shape and period of the DA pulses lead to the formation of an equilibrium, visible here as a plateau for the absolute amount of occupied receptor. This occurred at the level at which the amount of receptors unbinding until the next DA burst was the same as the amount of receptors getting occupied by the DA burst. Finally, the burst-pause events did not lead to an accumulation of occupied receptors over time. In fact, the receptor occupation was the same for burst-pause and non-event, except during the short burst component of the burst-pause events (Fig. 4, note the overlapping green and orange curves). This extends the property of burst-pause signals as “false alarm” signals to a wide range of occupancy levels.
Incorporating the slow kinetics in the model is crucial for functional considerations of the DA system. Currently, following the affinity-based model, the amplitude of a DA signal (i.e., peak [DA]) is often considered as a key signal, e.g., in the context of reward magnitude or probability (Morris et al., 2004; Tobler et al., 2005; Hamid et al., 2016). However, as DA unbinds slowly (over tens of seconds; Fig. 1d) and the binding rate changes approximately linearly with [DA], the amount of receptor occupancy does not primarily depend on the amplitude of the [DA] signal.
Because of the linearity of the binding rate, the receptor occupation increases linearly with time and [DA]–[DA]baseline, whereas the unbinding is negligible as long as t ≪ t1/2. Therefore the integral of the [DA] time course should be a close approximation of the receptor occupation for signals that are shorter than the half-life time of the receptors. We confirmed this consideration by simulating a range of DA signals (burst, burst-pauses, and ramps) with different durations and amplitudes. For each DA signal we compared its area under the curve with the resulting peak change in the absolute receptor occupancy. For both D1R and D2R we found that the maximum receptor activation was proportional to the area under the curve of the [DA] signal, although independent of its specific time course (Fig. 5). The small deviation from the proportionality seen for large-area DA signals for the D2Rs was due to the decrease in the amount of free receptor as more and more receptors were bound. In this regime the assumption that the binding rate is linear with [DA] was slightly violated leading to the non-proportionality. In contrast, the relationship between the DA burst amplitude and resulting receptor occupation was fit well by a quadratic function (Fig. 5b), which reflects quadratic increases in the area under the curve for larger amplitudes.
The overall striking proportionality of the integral of the DA signal with receptor binding indicates that D1Rs and D2Rs act as slow integrators of the DA signal. Interestingly, this means that DA ramps, even with a relatively small amplitude (Fig. 2b), are an effective signal to occupy DA receptors. In contrast, for locally very high [DA] (e.g., at corticostriatal synapses during phasic DA cell activity; Grace et al., 2007) our model predicts that the high concentration gradient would only lead to a very short duration of this local DA peak and thereby make it less effective in occupying DA receptors.
To further test the generality of our findings, we examined our model responses systematically for a set of different DA time courses (Figs. 6, 7). Although the shape of the DA pulses strongly affected the time courses of the receptor activation, D1 and D2 receptor activation were virtually identical for a given pulse shape. For DA bursts with different amplitudes (Fig. 6a), higher amplitudes of the DA burst lead to stronger receptor activation. However, the relationship between burst amplitude and receptor occupation was not linear, but instead reflected the area under the curve of the DA pulse (see previous paragraph).
Importantly, despite “slow” kinetics, the onset of the increase in [D1–DA] and [D2–DA] is immediate and determined by the area under the curve of the DA pulse up to each time point. This means that the DA receptor occupation reflects an ongoing integration of the DA signal. Furthermore, this provides us with an intuition for how the DA receptor occupation develops during a particular DA signal. The increase of the receptor occupation is controlled by the DA pulse shape and is proportional to the area under the curve of the DA pulse up to each time point. So with realistic kinetics the receptor occupation will always reach its maximum toward the end of the DA pulse and its rise profile depends on the specific pulse shape.
For a fixed burst amplitude, we also determined the effect of different DA re-uptake rates to look at potential differences in DA signaling in dorsal and ventral striatum, with fast and slow re-uptake, respectively. This was done by changing the parameter Vmax (see Materials and Methods), which controlled the time the [DA] took to return to the baseline from the peak value (Fig. 6b). Although this had only a small visible effect on the input DA signal (Fig. 6b, top), the resulting [D1–DA] and [D2–DA] were quite different. The difference in the resulting occupancy levels can again be understood in terms of the role of the pulse shape discussed above. For fast uptake, i.e., high Vmax, DA is cleared faster from the system, which reduces the area under the curve of the DA pulse but does not affect the peak DA levels during the burst. This leads to an important distinction between the affinity-based model and the model with realistic kinetics. With realistic kinetics the resulting receptor occupancy is lower for higher Vmax because the reduced area under the DA pulse gives the receptors less time to get occupied during high [DA] conditions (Eq. 5). By contrast, this property is not seen in the affinity-based model in which the time course of [D1–DA] and [D2–DA] simply follows the input [DA] signal, and thereby, the peak receptor occupation levels are not affected by Vmax.
Next, we examined DA ramps with different time courses, but the same maximal amplitude. Again, consistent with our consideration of the important role of the area under the curve of DA signals, we found that longer ramps lead to larger DA receptor occupation (Fig. 7a). We then investigated the DA signals that included the effects of pauses in DA cell activity further. First, we tested burst-pause signals and determined the role of the duration of the pause, when the [DA] decayed to zero. For a fixed burst amplitude and duration, a different duration of the subsequent pause lead to differing receptor activation levels when the burst-pause signal was over (Fig. 7b). This indicates that DA pauses are very effective in driving the receptor occupation quickly back to baseline (i.e., within few seconds) because, in this case, the receptor occupation changes reflect solely the unbinding rates. In contrast, for a burst followed by a return to baseline [DA], the decrease in receptor occupation would be slower because during the baseline portion of the signal both binding and unbinding processes play a role, and the binding counteracts some of the unbinding (Eqs. 5, 6).
In this context we also looked at pure DA pauses (i.e., without a preceding burst), e.g., reflecting DA cell responses to aversive stimuli (Schultz, 2007) that lead to reductions in [DA] (Roitman et al., 2008). These signals also lead to fast decreases in [D1–DA] and [D2–DA], with the duration of the pause having a strong effect on the amplitude and duration of the decrease (Fig. 7c). Pauses in DA cell firing may not necessarily lead to a [DA] of zero, as e.g., local mechanisms of DA release at terminals may persist even when DA neurons do not spike. Therefore, we also examined a scenario in which the DA pauses involved a reduction of [DA] to a quarter of its baseline level (Owesson-White et al., 2012). We found that in this case too, the pauses were effective in driving the receptor occupation back to the baseline (Fig. 7c, dashed line). However, the net unbinding was slower than during conditions when [DA] was zero during the pause. This is consistent with our corresponding observations for the burst-pause signals above (Fig. 3).
D1R and D2R occupancy in a probabilistic reward task paradigm
A general effect of the slow kinetics was that DA receptors remained occupied long after the DA pulse was over (Fig. 2), so that the effect of DA pulses was integrated over time (Fig. 4). To investigate the information that is preserved in the receptor occupation about DA signals on time scales relevant for behavioral tasks, we simulated sequences with probabilistic DA events (see Materials and Methods). First, we compared sequences, in which every 15 ± 5 s there was a DA burst with either 30, 50, or 70% probability (Fig. 8a,c). The resulting changes in [D1–DA] and [D2–DA] confirmed the integration of DA pulses over minutes. The integration of DA bursts was due to DA bursts arriving before the receptor occupation caused by the previous pulses had decayed, leading to an increased receptor activation compared with single DA bursts (Fig. 4). We then examined whether the DA receptor occupancy can distinguish different reward probabilities by using a simple classifier comparing two sequences with each other (see Materials and Methods). We tested sequences from 0 to 100% probability in steps of 10%, and ordered the resulting classification success in terms of the difference in reward probability between the two sequences (Fig. 8b,d,e,f). For example, a comparison between a 30 and 70% reward probability sequence yields a data point for a 40% difference. For both D1 and D2 receptors, we found that already for differences of 10% the classification exceeded chance level, and yielded near perfect classification an ∼40% difference. Overall, the classification was slightly better for D1R because of their slightly less developed plateauing response (Fig. 4). The successful classification of reward probabilities demonstrates that it would be possible for striatal neurons to read out different reward rates from DA receptor occupancy in a behavioral task. The classification in this example is only possible because of the slow kinetics of the receptors and, importantly, not due to an accumulation of DA. In this simulation there was no accumulation of DA itself and the [DA] dropped back to baseline in between the DA events. This provides a potential neural substrate for how fast DA signals can lead to an encoding of the slower reward rate, which can be used as a motivational signal (Mohebi et al., 2019).
Validation for fast binding kinetics
Our model assumption of slow kinetics was based on neurochemical estimates of wild-type DA receptors (Burt et al., 1976; Sano et al., 1979; Maeno, 1982). In contrast, recently developed genetically-modified DA receptors, used to probe [DA] changes, have fast kinetics (Patriarchi et al., 2018; Sun et al., 2018). Although the kinetics of the genetically modified DA receptors are unlikely to reflect the kinetics of the wild-type receptors (see Discussion), we also examined the effect of faster DA kinetics in our model. Fast kinetics were implemented by multiplying kon and koff by a factor q, keeping KD constant. We found that the similarity between [D1–DA] and [D2–DA] persists even if the actual kinetics were a 100 times faster than assumed in our model (Fig. 9). This was the case for all types of [DA] signals because the difference between the aggregate D1 and D2 binding rates (Eq. 5) still only differed by a factor of 1.5. Furthermore, the D2Rs did not show visible saturation effects even for q = 100. Faster kinetics mostly affected the amplitude of the receptor response and the time it took to return to baseline receptor occupancy. However, only for q = 100 the receptor occupation dropped slightly below baseline during the pauses of a burst-pause DA signal (Fig. 9a,b). On a longer time scale with repetitive DA bursts (Fig. 9e,f) D1Rs and D2Rs integrated the DA bursts over time even if kinetics were twice as fast (q = 2). This is because the half-time of the receptors were 40 s (for q = 2), while the DA burst signal was repeated every 15 s. Thereby, [D1–DA] and [D2–DA] were dominated by the repetition of the signal rather than by the impact of individual DA burst signals. In contrast, for q = 10 the change in receptor occupancy was dominated by the single pulses, because the half-life time was 8 s, whereby the receptors mostly unbind in between DA pulses. Therefore, our results concerning the similarity of D1 and D2 receptors do not depend on the exact kinetics parameters or potential temperature effects, as long as the parameter changes are approximately similar for D1 and D2 receptors. However, DA receptor kinetics faster by a factor of 10 or more affected the ability of DA receptor occupancy to integrate DA pulses over time (Fig. 9e,f).
In our model we assumed homogeneous receptor populations, namely that all D1 receptors have a low affinity and that all D2 receptors have a high affinity. However, this could be a simplification, as ≈10% of D2 receptors may exist in a low affinity state, while ≈10% of the D1 receptors may be in a high affinity state (Richfield et al., 1989a). Therefore, we also incorporated different affinity states of the D1 and D2 receptors into our model. The D1Rs in a high affinity state (D1high) were modeled by increasing the on-rate of the D1R but keeping its off-rate constant, creating a receptor identical to the D2high receptor. Although the high affinity state kinetics of the D1R are currently unknown, we choose this model as a faster on-rate potentially has the strongest effect on our conclusions. Correspondingly, we modeled the D2low receptor as a D2R with slower on-rate, which was equivalent to simply reducing [D2tot] because the D2low receptors were predominantly unoccupied during baseline DA and bound only sluggishly to DA during phasic signals. The main effect of incorporating the different receptor affinity states was a change in the respective equilibrium values of absolute concentration of receptors bound to DA (Fig. 10). However, importantly, taking into account these different affinity states, preserved the similarity of time courses of D1R and D2R occupancy and the ability to integrate DA pulses over time (Fig. 10, Fig. 11), because the half-life time of both receptors remained long.
Discussion
The functional roles of DA in reward-related learning and motivation have typically been studied by characterizing the firing patterns of dopaminergic neurons and the resulting changes in striatal [DA] (Schultz, 2007). In contrast to other, more conventional neurotransmitters like glutamate or GABA, the release of DA in the striatum may form a global signal that affects large parts of the striatum similarly (Schultz, 1998). Such global [DA] changes involve longer time scales lasting at least several seconds (Roitman et al., 2008; Howe et al., 2013). Importantly, to affect neural activity in the striatum, DA first needs to bind to DA receptors. This process is often simplified by assuming that this happens instantaneously, so that every change in [DA] is immediately translated into a change in DA receptor occupation. As this contradicts physiological measurements of the receptor kinetics (Burt et al., 1976; Sano et al., 1979; Nishikori et al., 1980; Maeno, 1982), we developed and investigated a model incorporating DA receptor kinetics as well as differences in D1 and D2 receptor abundance in the striatum.
Our results cast doubt on several long-held views on DA signaling. A common view is that D1 and D2 MSNs in the striatum respond to different DA signals due to the affinity of their predominant receptor type. Accordingly, phasic DA changes should primarily affect D1 MSNs, while slower changes or DA pauses should primarily affect D2 MSNs. In contrast, our model indicates that both D1R and D2R systems can detect [DA] changes, independent of their timescale, equally well. That is, slow tonic changes in [DA], phasic responses to rewards, and ramping increases in [DA] over several seconds lead to a similar time course in the response of D1 and D2 receptor occupation in our model. However, the baseline level of activated DA receptors and the amplitude of the response were typically twice as high in D2 compared with D1 receptors. Although, D1 and D2 receptors have opposing effects on the excitability (Flores-Barrera et al., 2011) and strength of cortico-striatal synapses (Centonze et al., 2001), we challenge the view that differences in receptor affinity introduce additional asymmetries in D1 and D2 signaling. Instead of listening to different components of the DA signal, D1 and D2 MSNs may respond to the same DA input. This would actually increase the differential effect on firing rate responses of D1 and D2 MSNs because the opposite intracellular effects of D1 and D2 activation (Surmeier et al., 2007) occur then for the whole range of DA signals.
Recently, ramps in [DA], increasing over several seconds toward a goal, have been connected to a functional role of DA in motivation (Howe et al., 2013; Hamid et al., 2016). In our model DA ramps were very effective in occupying DA receptors due to their long duration. In contrast, for brief phasic increases, the receptor occupation was less pronounced. Overall, our model predicts that the area under the curve of DA signals determines the receptor activation, which puts more emphasis on the duration of the signals, rather than the amplitude of brief DA pulses.
Our model is also relevant for the interpretation of burst-pause firing patterns in DA neurons. These are a different firing pattern than the typical reward-related bursts, and consist of a brief burst followed by a brief pause in action potentials. Such two-component responses of DA cells may reflect saliency and value components, respectively (Schultz, 2016). For example, during extinction learning burst-pause firing patterns have been observed as a response to conditioned stimuli, with each component lasting ∼100 ms (Pan et al., 2008). Our model provides a mechanistic account for how the burst-pause DA signals have a different effect on MSNs than pure burst signals, which is important to distinguish potential rewarding signals from other salient, or even aversive, stimuli. In our model the pause following the burst was very effective in reducing the number of occupied receptors quickly, thereby preventing the otherwise long-lasting receptor occupation due to the burst. Thereby, canceling the effect of the brief burst might be a neural mechanism to correct a premature burst response that was entirely based on saliency rather than stimulus value (Schultz, 2016). Because fast responses of DA cells to potentially rewarding stimuli are advantageous to quickly redirect behavior, the subsequent pause signal might constitute an effective correction mechanism labeling the burst as a false alarm.
Functionally, the slow unbinding rate of D1 and D2 receptors pointed to an interesting property in integrating phasic DA events over time. The unbinding rate might be one of the mechanisms translating fast DA signals into a slower time scale, which could be a key mechanism to generate motivational signals (Mohebi et al., 2019). Importantly, the slow kinetics of receptor binding do not prevent a fast neuronal response to DA signals. In our model [DA] changes affected the number of occupied receptors immediately; it just took seconds or even minutes until the new equilibrium was reached. However, reaching the new equilibrium is not necessarily relevant on a behavioral level. Instead the intracellular mechanisms that react to the receptor activation need to be considered to determine which amount of receptor activation is required to affect neural activity. In our model changes on a nanomolar scale occurred within 100 ms, a similar timescale as behavioral effects of optogenetic DA manipulations (Hamid et al., 2016).
The slower time scales were introduced into our model by the kinetics based on in vitro measurements (Burt et al., 1976; Sano et al., 1979; Nishikori et al., 1980; Maeno, 1982). A limitation of our model is the uncertainty about the accuracy of these measurements, and whether they reflect in vivo conditions. We addressed this here by also examining faster kinetics, for which there is currently no direct evidence in the literature. However, recently DA receptors have been genetically modified to serve as sensors for fast [DA] changes (Patriarchi et al., 2018), which suggests possible fast kinetics. It seems unlikely though that the kinetics of the genetically-modified receptors represent the kinetics of the wild-type DA receptors, as e.g., the screening procedure to find suitable receptor variants yielded a large range of different affinities (meaning changes in the kinetics of binding, unbinding or both) based on changes at the IL-3 site (Patriarchi et al., 2018). Changes in the IL-3 site have also previously been shown to strongly affect the receptor affinity (Robinson et al., 1994).
In addition to the receptor kinetics, the different abundances of D1 and D2 receptors are also key parameters in our model, which we estimated based on previous experimental studies. However, in case our estimates of the receptor abundances were incorrect, the receptor occupation would still be determined by the kinetic parameters and it would differ substantially from instant kinetics assumed in the affinity-based model. In particular, both receptor types would still not saturate during DA pulses, but integrate the DA signal over longer time scales. Furthermore, our results do not depend on the exact absolute receptor abundances, but on the relative abundances of D1 and D2 receptors. Therefore, our results on the similarity of the D1 and D2 responses hold as long as the abundance of the D1 receptors is approximately an order of magnitude higher than the abundance of the D2 receptors. Overall, we conclude that it is important to consider the effect of the receptor kinetics on DA signaling, which have not received much attention in experimental studies, nor in theoretical considerations of DA function thus far.
Footnotes
This work was supported by the University of Sheffield and its high performance computing resources, by funding from the EU H2020 Program as part of the Human Brain Project (HBP-SGA1, 720270; HBP-SGA2, 785907), the BrainLinks-BrainTools Cluster of Excellence funded by the German Research Foundation (Deutsche Forschungsgemeinschaft, Grant EXC 1086), and the state of Baden-Wuerttemberg through bwHPC. We thank Joshua Berke, Paul Overton, Alejandro Jimenez, Mohammadreza Mohagheghi Nejad, and Amin Mirzaei for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Lars Hunger at pc4xlh{at}sheffield.ac.uk