The dopamine dynamics model describes the temporal dopamine concentration, both tonic and phasic, autoreceptors occupancy and dopaminergic receptors occupancy. It was adapted from previously published models (Dreyer et al., 2010; Dreyer and Hounsgaard, 2013; Dreyer, 2014; Fuller et al., 2019).

The main mechanisms of dopamine regulation are outlined in the equations of the model and are represented in Figure 1. Dopamine is synthesized in the dopaminergic neurons and then released in the synaptic cleft. Sustained dopamine release refers to tonic dopamine, while transient dopamine release generated by bursts refers to phasic dopamine. The release of phasic dopamine is a reward prediction error signal (RPE) (Waelti et al., 2001; Marinelli and McCutcheon, 2014), whereas a drop in dopamine levels is a punishment prediction error signal. In the synaptic cleft, dopamine can be recaptured by DATs into the presynaptic neuron or be removed from the synaptic cleft by different mechanisms such as diffusion or inactivation by the Catechol-O-methyltransferase.

The remaining dopamine molecules can bind to dopaminergic autoreceptors located on the presynaptic neurons or to receptors on the postsynaptic neurons. In the present work, only dopaminergic receptors D₁ and D₂ are considered. All the above mentioned mechanisms are accounted for by the dopamine dynamics model, formulated in Equations (1) and (2), where C_DA(t) is the dopamine concentration (μM/L) in the synaptic cleft and AR(t) the autoreceptors occupancy.

As indicated in Equation (1), the release of dopamine is divided into two terms to account for both tonic and phasic release. The recapture by DATs is a saturable process described by a Michaelis-Menten equation. All other mechanisms contributing to dopamine removal are assumed to be linear (Budygin et al., 2002; Dreyer, 2014) and are schematized through the last term in the right-hand member of Equation (1). The binding to autoreceptors is proportional to dopamine concentration and free autoreceptors, while unbinding is proportional only to bound autoreceptors.

Autoreceptors have a regulatory effect on dopamine concentration. Indeed, they provide a negative feedback to adjust dopamine concentration through firing rate, synthesis, and release (Benoit-Marand et al., 2001; Beaulieu and Gainetdinov, 2011). Prolonged dopamine agonist exposure desensitizes autoreceptors in dopamine neurons (Robinson et al., 2017). Loss of inhibition influence facilitates further dopamine release and has been linked to drug abuse. Desensitization was not included in the model which is focused on the short-term effect of dopamine on autoreceptors. If tonic dopamine level decreases (in our ADHD model through increased dopamine reuptake), the temporary decrease in autoreceptor-mediated inhibition would mainly increase phasic dopamine release following the model developed by Grace (Grace, 1991, 2016). Autoinhibition of the presynaptic neurons is included in the model through the phasic release term only which is associated with the reward prediction error, while the tonic term is not here modified by autoreceptors occupancy (Grace, 1991).

The tonic dopamine release term is given by:

where ρ is the terminal density, Prtonic the tonic release probability, n₀ the number of molecules released per vesicles fusion, α_vf the extracellular volume fraction, N_A the Avogadros constant and υ_tonic the tonic firing rate. The tonic release is independent of autoreceptors occupancy, as explained above.

The phasic release term at time t is given by:

The terminal density (ρ), the number of molecules released per vesicles fusion (n₀), the extracellular volume fraction (α_vf) and the Avogadros constant (N_A) parameters are not modified by autoreceptors occupancy. Since vesicular release probability (Prphasic) and phasic firing rate (υ_phasic) are decreased by autoreceptors (Grace, 1991), they are assumed to be inversely proportional to autoreceptors occupancy (Beaulieu and Gainetdinov, 2011; Dreyer and Hounsgaard, 2013). The exact relationship is not known but assumed here as inversely proportional for simplicity. The value 0.334, used to normalize the equation for the control case, corresponds to autoreceptors occupancy. Therefore, Equation (5) indicates that the activation of autoreceptors reduces phasic dopamine release. The values 0.1s (Bamford et al., 2018) and 0.05s represent the latency and duration of the reward or punishment error prediction signal, respectively. Phasic dopamine release is also proportional to the reward prediction signal (RPE). This issue will be discussed in more details in Section 2.3.

In the occurrence of a punishment, the activity of the dopamine neuron is temporarily suppressed (both tonic and phasic firing rate fall to zero). According to Equations (1) and (3), this can be simulated in the model assuming υ_tonic = 0 which corresponds to the following differential equation:

With the parameters we used, this equation requires about 500 ms to reach the new equilibrium with C_DA = 0, which is close to the duration of dopamine neuron activity suppression after the absence of an expected reward (Schultz et al., 1997). However, the time to reach this equilibrium may vary as a function of the previous discharge rate, tonic dopamine level, or reuptake. To simplify the model, the value C_DA = 0 was directly applied at the same time as for the phasic dopamine discharge associated with a reward, as shown in Equations (7) and (8). Setting the dopamine concentration at zero instantaneously when a punishment occurs is a simplification of the physiologic mechanisms and the pause in the firing rate was defined as in Dreyer et al. (2010). This simplification was used since the purpose of this work was to study the behaviors in a qualitative manner. In future work, we will implement more physiologic parameters with their variability.

In the model, autoreceptors occupancy depends on the overall dopamine concentration (tonic and phasic). It could be argued that, due to diffusion, only a fraction of phasic dopamine reaches autoreceptors and thus alters the release. Simulations were performed to integrate this concentration gradient on phasic dopamine reaching autoreceptors, but the results were not significantly different (not shown here), therefore the version presented here was chosen for simplicity.

Finally, dopamine molecules can bind to dopaminergic receptors, corresponding to D₁ and D₂ receptors in the current work. The occupancy of receptors of type i ∈ {1, 2} in time is given by the following equation:

where BmaxDi and kDDi are the maximal concentration and dissociation constant of type i receptors, respectively. Receptors occupancy will be used in the neurocomputational model of basal ganglia as the postsynaptic effect of dopamine on the neurons in the different neurotransmission pathways (Hille, 1992).

The parameter values for the dopamine dynamics model are given in Table 1. As mentioned in this Table, the dopaminergic terminal density was adapted. As this density is inhomogeneous (Dreyer, 2014; Fuller et al., 2019), its value was set to obtain a tonic dopamine concentration in the control group of 0.02 μM/L as reported in the literature (Wanat et al., 2009; Hunger et al., 2020).

Using the developed model, two groups of virtual individuals were created: control and dopamine imbalance individuals. The difference between the two groups lies in the modification of the V_max parameter of Equation (1). From a mathematical standpoint, the parameter k_m could also have been decreased to obtain similar results.

Tonic and phasic dopamine are coding prediction error signals in the basal ganglia (Schultz, 2017). ADHD is associated with dopamine dysfunctions in the cortex and the basal ganglia (Giedd et al., 2001; Seidman et al., 2005; Nakao et al., 2011; Cubillo et al., 2012; Frodl and Skokauskas, 2012; Oldehinkel et al., 2016). Hence, a neurocomputational model of basal ganglia with a learning procedure was added to the dopamine dynamics model.

The neurocomputational model presented here is an adaptation from the model developed in Baston et al. (2016). It involves the temporal neural activity in the cortex, the thalamus and the different regions of the basal ganglia (striatum, globus pallidus pars interna and pars externa, and subthalamic nucleus), with a representation of the external stimulus S. The neuronal activities are normalized to obtain a value between 0 and 1. The connection between each region follows three neurotransmission pathways: direct, indirect and hyperdirect. The direct pathway promotes movement, the indirect inhibits it, and the hyperdirect pathway suppresses erroneous movements. D₁ and D₂ receptors occupancy have an excitatory effect in the direct pathway and an inhibitory effect in the indirect pathway, respectively. Both pathways are potentiated by the effect of cholinergic interneurons, also included in the model.

A representation of the neurocomputational model of basal ganglia is given in Figure 2. Each region of the model is divided into four action channels, representing different alternative choices. This division allow investigating the response of basal ganglia to various target stimuli. Neural activity in each action channel is computed through an ordinary differential equation, simulating neural dynamics, and a sigmoidal relationship, which mimics the typical non-linear phenomena of the neurons (lower threshold and upper saturation). The input to each differential equation is calculated by summing all the upstream activities converging to that neuron, weighted by the synaptic strength. The synaptic weight matrices correspond to the weight of connections between the regions for all four action channels. Equations and parameter values of the model can be found in the Supplementary Material.

Impairments in reinforcement learning are thought to be involved in ADHD (Sagvolden et al., 2005; Tripp and Wickens, 2008; Alexander and Farrelly, 2018). Therefore, we included a reinforcement learning process with reward and punishment prediction error signals in the model. The strength of connections between each region of basal ganglia is given by synaptic weight matrices noted w^ij, where i and j are the postsynaptic and presynaptic regions, respectively. The values of these weights can be modified by the learning process. For simplicity, only matrices related to striatum, w^GS, w^NS, w^GC, w^NC, were considered to be plastic; these connections are represented by dashed lines in Figure 2. The matrices w^GC and w^NC are diagonal while w^GS and w^NS are full matrices. At the beginning of the learning process, these weight matrices are in a naive state, with no differentiation between the actions channels. Here are the initial value of the matrices:

We here give the details of a typical trial of the learning process. A stimulus representation S is sent for 800 ms to each action channel. One channel will receive a strong stimulus of value 1, another one receives a weaker stimulus of value 0.2, while the two others receive even weaker stimuli with a value of 0.1 each. In the present work, we used an input vector with the same dimension as the number of possible actions, with a higher value (close to 1) at the same position of the rewarded action, and a smaller value at the positions of the punished actions, just to simplify the final analysis of the synapses. An input vector with different dimensions and with different values could be used as well, resulting in a more complex pattern of synapses. The idea here is to simply associate an input vector to a “winner takes all” output vector, considered as the selected response. The possible considered vectors for S are S=[10.20.10.1], S=[0.210.10.1], S=[0.10.10.21] and S=[0.10.110.2]. Neuronal activity in all regions of basal ganglia are computed for 800 ms. An action is considered to have been performed or chosen if the activity in its related action channel in the cortex is above 0.9, while the activity in all other channels is close to zero, using the winner-takes-all dynamics implemented in the cortex.

We used a fixed scale of prediction error values throughout learning. The prediction error is the discrepancy between observed and expected outcome, and a naive subject cannot predict whether the response would be correct or not. If the chosen action is in the action channel with the highest value of S, a reward prediction error of 1 is attributed. If however the second highest value (0.2) is chosen, a smaller reward prediction error of 0.1 is attributed. A punishment prediction error is given when the lowest value (0.1) is chosen. Rewards prediction errors are signaled by phasic dopamine peaks governed by Equation (6). When a punishment prediction error occurs, dopamine concentration drops to zero. This is equivalent to providing the virtual subjects with rewards and punishments, but we delivered directly the reward/punishment prediction error dopamine signals. These prediction errors are the differences between received and predicted rewards (Schultz, 2016), where here the virtual patient always predicts a reward when an action is chosen. This process is repeated over 1,000 trials (epochs). Once the learning procedure is complete, the model is expected to effectively differentiate between weak and strong stimuli, so that responses occur only when strong stimuli are applied.

The resulting rewards/punishments prediction error signal will lead to a modification of the synaptic weights contained in the matrices. These weights modifications during the learning process are dictated by the Hebb Rule, which states two neurons having both high activity will strengthen their connection, whereas connection will weaken in case of neurons with opposite activity. The Hebb rule describes how much the weights are increased or decreased at each step of the training procedure. In particular, the following equation holds at each step to assign a new synaptic value, Baston and Ursino (2015):

where w^AB represents the matrix containing the weights from the presynaptic region B to the postsynaptic region A, with A being either S or C in Figure 2 and B being either G (Go) or N (NoGo) in the same figure, and Δw^AB is the synaptic change computed at that step. Each row in these matrices represent the synapses entering the postsynaptic neuron, and each column those emerging from the presynaptic one. Hence, all matrices have 4 × 4 dimensions in the work presented here. This modification of the synaptic weights happens once every epoch between a latency period of 0.1s and for a duration of 0.05s once an action is chosen. The latency and duration are the same as the ones for the reward/punishment error prediction signal. The individual elements at position ij in the array Δw^AB are computed through the following equation (Hebb rule):

where yjB is the activity of the presynaptic neuron in the action channel j of the region B, yiA is the activity of the postsynaptic neuron in the action channel i of the region A and ϑ_presynaptic, ϑ_postsynaptic the pre- and postsynaptic thresholds. The positive part function ([]⁺) ensures that learning occurs only if the presynaptic neurons are excited and their activity is above the threshold. Dopamine is thought to have the ability to modulate synaptic plasticity, although the exact relationship does not seem to be established (Reynolds and Wickens, 2002; Frémaux and Gerstner, 2015; Madadi Asl et al., 2019). From previous work, it seemed reasonable to assume a proportional relationship with dopamine ratio and RPE. Of course, in case of diagonal matrices (w^GC and w^NC), only the elements with i = j are trained, compared to non-diagonal matrices w^GS and w^NS where all elements are trained. The gain parameter ϕ is proportional to the reward prediction error since, for example, a large reward prediction error will lead to a larger variation in the synaptic value than a small reward prediction error. The gain parameter is also proportional to the ratio of phasic peak and tonic dopamine. This ratio is calculated beforehand and considered as a constant. The equation is the following:

The dopamine ratio is higher in the dopamine imbalance group (with a value of ~ 8.3) compared to the control one (with a value of ~ 3), so the gain parameter ϕ is higher.

The control and dopamine imbalance groups, with 10 virtual subjects each, were created with the model. The only difference between the two groups is in the value of V_max. A higher rate of dopamine recapture is expected to lower the dopamine tonic concentration which in turn is expected to increase the phasic dopamine concentration, and thus in the tonic phasic dopamine ratio, through a lower occupancy of autoreceptor. The steps of the learning procedure of a subject are summarized below.

During the test phase, we also used a four-choice reaction time task. A series of stimuli are presented to the virtual individuals in the different action channels through a signal S of the neurocomputational model of BG to the cortex. The stimulus in the targeted action channel has a value of 1 with the addition of noise. Noise is also added in the other action channels directly in the cortex. Each stimulus is presented for 1, 800 ms with a 500 ms pause in between each stimulus. The criterion for a response is an activity in one of the four action channels in the cortex C, which constitute the output of the model, greater than 0.9. Due to the winner-takes-all dynamics, the other three channels will then have activity close to zero. For simplicity purposes, a response in the same action channel as the target stimulus is considered as a success. Otherwise it constitutes a failure. Or course successes and failures could have been defined in different ways. The idea here is simply to associate to an input vector, an output vector considered as the correct responses.

During the test phase, there is always a response after a stimulus, being a success or a failure. The number of correct answers or successes represents the performance of the virtual individuals. Each individual is presented 100 stimuli. The mean and standard deviation of the percentage of successes and of the reaction times are computed in each simulated group. Stimulus of different amplitudes were also sent in the first action channel and the responses were recorded to study the differentiation between weak and strong signals. In order to compare the ability of differentiating between weak and strong signals, we repeated the task and computed the cortex activity for different values of noise added to the input signal (S).

During the test phase, reaction times were also computed. The reaction time is here defined by the difference between the time at which the neuronal activity in one of the action channels reaches a value of 0.9 and the time at which the stimulus was sent in the sensory representation S.

Using the model, dopamine concentrations were simulated for the two groups as shown in Figure 3. Phasic peaks were created by a burst lasting 0.05 s.

As seen in Figure 3, dopamine imbalance individuals have lower tonic dopamine concentration due to higher dopamine recapture. In turn, autoreceptors regulation causes higher phasic dopamine concentration. This dopamine imbalance will have different impact on the learning process in the basal ganglia.

During the training phase, we computed the number of trials to obtain 5 successful responses over 10 successive trials. All participants in the normal group reached the learning criterion, but 2 participants in the dopamine imbalance group failed to do so even after 1,000 trials. The number of trials to reach criterion was on average 65.1 (SD = 52.6) in the control group, but 20% higher in the dopamine imbalance group, with an average of 85.5 (SD = 67.8), excluding those who never reached the criterion.

In the first task, the mean and standard deviation of the percentage of successes to a series of 100 stimuli and of reaction times are computed in each simulated group and shown in Figure 4.

The mean reaction time in the control group is 148 ms and the standard deviation is 5 ms. The mean percentage of successes is 100 with a standard deviation of 0. In the dopamine imbalance group, the mean reaction time is 166 ms with a standard deviation of 57 ms. The mean percentage of successes is 78 with a standard deviation of 19. As shown in Figure 4, the rate of successes was lower and more variable in the dopamine imbalance group, as compared to the control group. Moreover, the simulated mean reaction times was slower in the dopamine imbalance group than in the control group. In our simulations, the mean and standard deviation of reaction times are respectively, 1.12 and 11.4 times larger in the dopamine imbalance group than in the control group. The significance of the reaction time difference was not evaluated because only 10 patients were simulated in this study to present the model. Also, as described further, the patients in the dopamine imbalance group are heterogeneous and can be divided into three subgroups with different mean reaction times.

We used the ex-Gaussian distribution to estimate the reaction time distribution by combining a normal and an exponential distribution. Three parameters characterized the ex-Gaussian distribution: the mean μ and standard deviation σ of the normal distribution, and τ representing the mean and standard deviation of the exponential part. An ex-Gaussian distribution was fitted to the simulated reaction times of the virtual individuals as seen in Figure 5.

The τ parameter was 12 times larger in the dopamine imbalance group than in the control group (47 vs. 3.8) while the μ parameter was 0.82 times smaller (118 vs. 144).

We assessed the performance of the individuals in each group described in the above section by increasing the standard deviation of the noise added to the input signal S. A series of 100 stimuli was again presented with noise directly added to the stimulus representation in the cortex S, with a mean of 1 and a standard deviation ranging from 0 to 1. As the standard deviation of the noise increases, the probability of having high intensity noise increases which further complicates decision making for the virtual patients and therefore affects the percentage of successes. Figure 6 shows that in the dopamine imbalance group the mean percentage of successes (orange solid line) quickly dropped while the variability (orange shaded area) increased with increasing noise variability. By contrast, in the control group, the performance remained optimal, with no variability, until the noise variability was greater than 0.6.

During the test phase, we also computed the output activity in the cortex related to the response as a function of the input value of the stimulus. A stimulus of different amplitudes, ranging from 0.1 to 1, is sent in the first action channel while all three other channels receive noise of small amplitude. The mean, the 5th and the 95th output curves of the cortex neuronal activity in the first action channel as a function of the input signal value for each group are shown in Figure 7.

By comparing neural activity at basal ganglia input and output, it is clear that in control subjects, the basal ganglia have a high neural gain. Response-related activity is suppressed until stimulus-related cortical activity reaches 0.5 in the control group. Output activity then increases rapidly for an input between 0.5 and 0.7 at which point it remains maximal. In contrast, in the dopamine imbalance group, activity is suppressed up to an input of 0.4, after which the gain increases rapidly but only for stimulus-related activity between 0.4 and 0.5. For stimulus-related activity values between 0.5 and 1, the gain is strongly attenuated as response-related activity increases from 4.5 to 7. However, the most striking aspect of the gain is the extreme variability of the output in the dopamine imbalance group, which ranges from 0 to 1 in response to stimulus-related activity values between 0.7 and 1. In this group, some individuals respond correctly and others have wrong responses which will lead to an output activity close to zero due to the winner-takes-all dynamic, thus inducing high variability. In contrast, in the controls, the variability is almost zero, except for the amplification phase, especially around the inflection point.

Four synaptic weights matrices were modified during training: w^GS, w^NS (stimulus-related synaptic weights) and w^GC, w^NC (response-related synaptic weights). These matrices start in a naive configuration, with no differentiation between the four action channels. They are modified during the training by using the Hebb Rule, with a gain parameter that is proportional to the phasic vs. tonic dopamine ratio.

Over the course of the 1,000 trials in the training phase, the matrix weights changed differently between the two groups, and between individual subjects within each group. Indeed, the trends of synaptic weight evolution were the same for the control and dopamine imbalance groups, but inter-individual differences in synaptic weights and their evolution during learning were much larger in the dopamine imbalance group. Hence, inter-individual differences were much larger at the end of the learning phase in the dopamine imbalance than in the control group. More details on the evolution of the synaptic weight matrices are given in the Supplementary Material.

In the present section, a metric is developed to differentiate the performance in the test phase of the dopamine imbalance group from the control one based on their history during the training phase. During the training process, the history of rewards and punishments is stored in a vector with value 1 for a large reward, 0.1 for a small reward, −1 for a punishment and 0 for no response. It is therefore possible to study the history of each individual and to relate it to his performance in the test phase.

Figure 8 shows the cumulative sum of the history vector for each action channel of the first 5 individuals in each group. A negative cumulative sum results from a series of failures overcoming successes, while a positive cumulative sum would indicate the opposite.

There seems to be an initial phase in which there is an excess of errors. The virtual individuals start in a naive state, meaning no differentiation between the action channels. Hence, the initial responses have a random success rate of 25% and can lead to an excess of errors. In the second phase (> 500 epoch), rewards prediction errors dominate over punishments for all actions.

Individuals from the control group seem to learn each action in a proportional way for all action channels. The individuals in the dopamine imbalance group had a higher number of rewards for some action channels at the expense of the others. In order to quantify the inter-individual differences in learning, a weighted standard deviation (weighted std) for the cumulative sum of history was computed for each individual, and expressed by the following equations:

where i is the epoch number, j the action number, cumsum_actionj (i) the cumulative sum of history vector for action j at epoch i and mean(i) is the mean of cumulative history at epoch i for all action channels. The standard deviation of the history (std_history) is weighted by a ratio to take into account the fact that the cumulative sum of history is either positive or negative. The ratio is the sum of negative cumulative sum of history divided by the sum of positive cumulative sum of history, leading to a larger ratio when the negative cumulative sum exceeds the positive one. Division by 1,000 is for scaling. The weighted std_history was larger in the dopamine imbalance group than the control one. In order to assess the relationship between the training and test phase, a plot of the standard deviation of the reaction times as a function of the weighted std_history value is depicted in Figure 9. A linear regression (dashed line) and a quadratic function (dashed curve) between the weighted std_history and the standard deviation of reaction times were applied to the control group and the imbalance group, respectively. The individuals in the dopamine imbalance group could be divided into three subgroups (a, b, and c) along the quadratic regression as seen in Figure 9. Group a contained the individuals with a perfect performance, low μ, low σ and low τ, which explains their proximity to the individuals in the control group. The individuals less than perfect performance were divided into groups b (75% of successes) and c (60% of successes). The distribution of reaction times in the group b is closer to an exponential distribution than to a normal one with low μ and σ but very high τ. These individuals have both fast and very slow reaction times, driving thus the mean to a high value. As the weighted std_history increases for individuals in group c, the performance further decreased with fewer correct responses, the μ parameters increased, and the σ and τ had intermediate values and were quite similar.