Dopamine Neurons Change the Type of Excitability in Response to Stimuli

The dynamics of neuronal excitability determine the neuron’s response to stimuli, its synchronization and resonance properties and, ultimately, the computations it performs in the brain. We investigated the dynamical mechanisms underlying the excitability type of dopamine (DA) neurons, using a conductance-based biophysical model, and its regulation by intrinsic and synaptic currents. Calibrating the model to reproduce low frequency tonic firing results in N-methyl-D-aspartate (NMDA) excitation balanced by γ-Aminobutyric acid (GABA)-mediated inhibition and leads to type I excitable behavior characterized by a continuous decrease in firing frequency in response to hyperpolarizing currents. Furthermore, we analyzed how excitability type of the DA neuron model is influenced by changes in the intrinsic current composition. A subthreshold sodium current is necessary for a continuous frequency decrease during application of a negative current, and the low-frequency “balanced” state during simultaneous activation of NMDA and GABA receptors. Blocking this current switches the neuron to type II characterized by the abrupt onset of repetitive firing. Enhancing the anomalous rectifier Ih current also switches the excitability to type II. Key characteristics of synaptic conductances that may be observed in vivo also change the type of excitability: a depolarized γ-Aminobutyric acid receptor (GABAR) reversal potential or co-activation of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPARs) leads to an abrupt frequency drop to zero, which is typical for type II excitability. Coactivation of N-methyl-D-aspartate receptors (NMDARs) together with AMPARs and GABARs shifts the type I/II boundary toward more hyperpolarized GABAR reversal potentials. To better understand how altering each of the aforementioned currents leads to changes in excitability profile of DA neuron, we provide a thorough dynamical analysis. Collectively, these results imply that type I excitability in dopamine neurons might be important for low firing rates and fine-tuning basal dopamine levels, while switching excitability to type II during NMDAR and AMPAR activation may facilitate a transient increase in dopamine concentration, as type II neurons are more amenable to synchronization by mutual excitation.


Abstract
The dynamics of neuronal excitability determine the neuron's response to stimuli, its synchronization and resonance properties and, ultimately, the computations it performs in the brain. We investigated the dynamical mechanisms underlying the excitability type of dopamine (DA) neurons, using a conductance-based biophysical model, and its regulation by intrinsic and synaptic currents. Calibrating the model to reproduce low frequency tonic firing results in N-methyl-D-aspartate (NMDA) excitation balanced by γ-Aminobutyric acid (GABA)-mediated inhibition and leads to type I excitable behavior characterized by a continuous decrease in firing frequency in response to hyperpolarizing currents. Furthermore, we analyzed how excitability type of the DA neuron model is influenced by changes in the intrinsic current composition. A subthreshold sodium current is necessary for a continuous frequency decrease during application of a negative current, and the low-frequency "balanced" state during simultaneous activation of NMDA and GABA receptors. Blocking this current switches the neuron to type II characterized by the abrupt onset of repetitive firing. Enhancing the anomalous rectifier Ih current also switches the excitability to type II. Key characteristics of synaptic conductances that may be observed in vivo also change the type of excitability: a depolarized γ-Aminobutyric acid receptor (GABAR) reversal potential or coactivation of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPARs) leads to an abrupt frequency drop to zero, which is typical for type II excitability. Coactivation of N-methyl-D-aspartate receptors (NMDARs) together with AMPARs and GABARs shifts the type I/II boundary toward more hyperpolarized GABAR reversal potentials. To better understand how altering each of the aforementioned currents leads to changes in excitability profile of DA neuron, we provide a thorough dynamical analysis. Collectively, these results imply that type I excitability in dopamine neurons might be important for low firing rates and fine-tuning basal dopamine levels, while switching excitability to type II during NMDAR and AMPAR activation may facilitate a transient increase in dopamine concentration, as type II neurons are more amenable to synchronization by mutual excitation.

Introduction
Midbrain dopamine (DA) neurons predominantly fire in a low frequency, metronomic manner (i.e. tonic) and display occasional, yet functionally important, high frequency, burst-like episodes [1,2]. While regular tonic firing is observed in isolated preparations (i.e. slices), tonic firing pattern in vivo is somewhat more variable due to active synaptic inputs [3,4]. Such tonic activity is important for maintaining a constant basal level of dopamine in projection areas. Accordingly, abnormal basal DA levels are linked to psychiatric disorders from depression to schizophrenia [5,6]. While the maintenance of basal DA levels seem to be critical for normal brain function, a consistent picture has not yet emerged regarding how changes in firing patterns of the DA neuron facilitates this important biological function.
Background activity of the DA neuron appears to rely on the intrinsic pacemaking mechanism that generates tonic firing. The current composition producing low-frequency pacemaking in DA neurons is of vibrant debate among researchers in the field. A number of experimental [7][8][9][10][11][12][13][14][15][16][17][18][19] studies suggests that the maintenance of tonic firing in at least a subpopulation of DA neurons relies on the interactions of the voltage gated calcium (Ca 2+ ) and SK-type Ca 2+ -dependent potassium (K + ) currents Slow pacemaking in our model relies on a subthreshold Ca 2+ -K + oscillatory mechanism, similar to a number of well-established models [4,[20][21][22][23][24]. Interaction between Ca 2+ and Ca 2+ -dependent K + currents periodically brings the neuron to the spike threshold and generates a metronomic firing pattern. In our model, spike-producing currents (fast sodium and the delayed rectifier potassium) play a mostly subordinate role in this dynamic, adding a spike on top of the oscillations without significant changes to the period or shape of voltage and calcium oscillations, as in the study by Wilson and Callaway 2000 [4]. A number of studies suggest an additional Ca 2+ -independent oscillatory mechanism [25][26][27]. In particular, they emphasize the contribution of sodium currents to pacemaking. We review the literature on the mechanisms of DA neuron pacemaking in the discussion section in more detail. The specific composition of currents contributing to oscillations determines the response of the DA neurons to stimuli, their synchronization properties and, ultimately, the computations they perform. In this paper, we use recent experiments to calibrate the dynamical properties of the DA neuron and determine its excitability type.
A standard method to classify neuronal excitability is via characterizing the frequency-toinput relationship, or F-I curve. Two major types of excitability can be determined based on how the onset of tonic firing occurs as the applied current increases and the neuron is released from quiescence at the hyperpolarized rest state [28]. A type I-excitable neuron can fire at an arbitrary low frequency near the onset of firing, whereas a type II neuron shows a discontinuous jump to a minimal frequency above a certain current threshold and fires only in a limited range of frequencies [21]. Generally, the onset of repetitive firing occurs through one of two mechanisms: 1) a saddle-node bifurcation on invariant circle or SNIC (type I excitability) or 2) an Andronov-Hopf bifurcation (type II excitability) [29]. Type II neurons, such as fast-spiking inhibitory interneurons in the cortex, display precise spike timing even in the presence of noise and are therefore suitable for the implementation of spike time coding [30,31]. A type I neuron, such as a weakly adapting cortical pyramidal neuron, was shown to relay the stimulus rate by modulating its own frequency, and, therefore, displayed rate coding [31]. Further, type II neurons display resonance and controlled synchronization in networks [31][32][33]. For neurons that are tonically active without any injected current, such as DA neurons, the transition to the non-spiking rest state occurs when a sufficiently strong hyperpolarizing current is injected. For these neurons, the excitability type would be defined by the transition from tonically firing to quiescent/excitable: again type I would show a smooth frequency decrease to zero, while type II should show an abrupt transition to quiescence. At this point, there is no direct evidence defining to what type of excitability DA neurons belong, and how the different intrinsic conductances and the different synaptic inputs influence their type. Determining the type of excitability will allow us to predict the behavior of the DA neuron during application/blockade of different currents and better understand computations it performs in different input conditions (e.g. rate coding vs. resonance at a particular input frequency).
A number of experimental studies provide indirect evidence of excitability type of the neuron in control conditions and during activation of synaptic inputs. It has been shown that the firing rate of DA neurons increases linearly in response to a ramping depolarizing current until it goes into depolarization block (e.g. [3,34]). Further, injection of a tonic hyperpolarizing current to the regularly firing DA neuron in vitro increases its interspike intervals [7]. The firing properties of the neuron in response to a combination of tonic inhibitory and excitatory synaptic conductances were investigated by Lobb and colleagues [35,36]. Using the dynamic clamp technique, they injected inhibitory γ-Aminobutyric acid (GABA) and excitatory Nmethyl-D-aspartate (NMDA) receptor conductances in SNc DA neurons. Injection of tonic GABAR conductance decreased the firing rate of the neuron several-fold. Furthermore, the neuron fired at low frequencies when NMDAR and GABAR conductances balanced each other. Thus, NMDAR activation, which strongly increases the firing frequency [37][38][39][40], can be effectively compensated by GABAR activation. Such compensation would be impossible if the inhibition produced an abrupt transition to quiescence and the neuron jumped from a high frequency to zero. The experimentally observed compensation suggests, again, a smooth frequency decrease upon GABAR activation rather than an abrupt transition to the resting state at hyperpolarized potentials. Together, these data resemble the tonic firing/quiescence transition in type I neurons with two distinctions. First, the transition parameter is not an injected current, but an ohmic GABAR conductance. In experiments, a conductance has already been used instead of an injected current to determine the neuronal excitability [41]. Second, the coactivation of the NMDA receptor introduces an additional parameter (its maximal conductance). Both of these extend the definition of excitability into the space of synaptic conductances. Formally, the excitability type is an intrinsic property of a neuron, yet viewing synaptic inputs as changing excitability of a neuron is a powerful concept used to understand neuron dynamics in vivo [31,42]. We used the experiments described above to parameterize a model of the DA neuron, determine its type of excitability, and determine how intrinsic and synaptic currents shape the excitability type and, therefore, the computational properties of the neuron.
These experiments suggest that the DA neuron exhibits type I excitability in isolation from synaptic inputs and under the balanced influence of excitatory and inhibitory synaptic conductances. However, the excitability type has been shown to vary depending on the intrinsic currents and network connectivity [42,43]. For example, modeling results suggest that changes in the intrinsic currents, e.g. L-type Ca 2+ current, can switch the excitability type of the DA neuron [44,45]. Here we address the variability in the excitability type under different conditions by studying the contribution of intrinsic and synaptic currents to regulation of the lowfrequency DA neuron firing.

Compensatory action of asynchronous NMDA and GABA inputs
We investigated the behavior of a simulated dopamine neuron in response to irregular asynchronous GABA and glutamate (Glu) inputs to mimic temporal structure of neural firing in in vivo conditions. The Glu input was produced by Poisson distributed spike trains and GABA inputs was explicitly modeled as activity of a population of GABA neurons (detailed description of the inputs and equations are given in the methods section). We quantified changes in the firing rate and the regularity of DA neuron firing in response to synaptic inputs of different strengths (Fig 1A1 and 1A2). We identified a parameter region where the excitatory and inhibitory inputs balance to produce low frequency DA neuron firing at rates similar to background firing (Fig 1A1, between the black lines). This happens because asynchronous GABA and Glu inputs (see rasters in Fig 1B1 and 1B2) activate GABARs and NMDARs nearly tonically ( Fig  1B3 and 1B4) and provide quasi-constant levels of inhibition and excitation to the DA neuron respectively. Under the influence of these two inputs, the DA neuron fires similarly to the in vitro-like conditions (tonic inputs), but with less regularity, which is typical of the background firing in vivo. An example voltage trace of the DA neuron in response to synaptic inputs formed by asynchronous Glu and GABA populations is shown in Fig 1B5. Fluctuations in the firing of neural populations innervating the DA neuron can produce irregular spiking as observed in in vivo experiments.
Considering that asynchronous inputs produce nearly constant receptor activation (see Fig  1B3 and 1B4), for further analysis we substituted these Glu and GABA inputs by tonic currents. Moreover, tonic synaptic currents mimic long-lasting injection of the conductances in dynamic clamp experiments [46,35,36], or iontophoresis of the agonists [39,40], or bath application of the agonists. Transition from asynchronous inputs to tonic currents is described in the methods section.

Balance of tonic NMDA and GABA inputs
Our next goal was to reproduce the experimentally-observed compensatory influence of tonic NMDAR and GABAR conductances [35]. Using the dynamic clamp technique, it was shown in vitro that a balanced injection of GABAR and NMDAR conductances leads to DA neurons firing at frequencies comparable with background frequencies (1)(2)(3)(4)(5). Removal of inhibition in such conditions evokes a classical disinhibition burst (the disinhibition model of burst generation is well known and described in e.g. [35,[47][48][49] [35]. In this example, the simulated DA neuron is tonically active at 1.5 Hz during tonic co-activation of NMDA and GABA receptors (g NMDA = 16.9 m S/cm 2 , g GABA = 5mS/cm 2 ). Removal of the GABAR conductance produces an episode of high-frequency firing (Fig 2A). Removal of the NMDAR conductance produces a pause in firing (Fig 2A).
We explored the range of NMDAR and GABAR conductances that produce tonic firing in the DA neuron model (Fig 2B). Compensation of NMDAR and GABAR activation can be readily achieved near the upper boundary of the firing region (Fig 2B, blue). When both receptors are activated, low frequency tonic activity is observed (Fig 2A). The dot labeled as A on the heat plot indicates conductances taken for this simulation. As in the experiments [35], the balanced region is stretched linearly on the conductance plane with NMDA/GABA slope around 3.4. Moving to the left on the diagram corresponds to deactivation of the NMDAR current and blocks DA neuron firing due to the remaining GABAR activation (Fig 2A). A pause may also be produced by stronger activation of the GABAR (Fig 2C). Conversely, moving down on the diagram corresponds to deactivation of the GABAR and evokes high-frequency firing (Fig 2A). The firing frequency also increases by moving from the upper boundary of the firing region to the right (increasing NMDAR conductance; Fig 2D). These two directions correspond to two ways of eliciting a DA neuron burst: strengthening NMDA excitation or removing inhibition, respectively.
Excessive tonic NMDAR activation leads to a depolarization block, as shown in Fig 2B and  2E at high NMDA and low GABA receptor conductances. Interestingly, application of a tonic GABAR conductance in combination with an excessive NMDAR conductance may rescue high-frequency firing in the model (Fig 2E). Thus, the compensatory influence of GABAR activation removes depolarization block induced by an excessive NMDAR activation and restores the intrinsic oscillatory mechanism required for tonic firing. The smooth frequency decrease to zero as the neuron transition to quiescence when GABAR conductance increases suggests type I excitability for the DA neuron both in in vivo and in vitro like cases (Figs 1A1 and 2B). However, excitability is classically defined by the structure of the transition between spiking and hyperpolarized rest state induced by an injected current, as opposed to a synaptic conductance. We show that the DA neuron exhibits type I excitability by standard definition with a continuous F-I curve and place it in Supporting Information (S1 Fig) as this case has less physiological significance than the influence of synaptic currents. We further investigate the influence of intrinsic and synaptic currents on the excitability type of the DA neuron.

Influence of intrinsic currents on the type of excitability
The role of Ca 2+ and Ca 2+ -dependent K + currents. The subthreshold Ca 2+ -K + oscillatory mechanism underlies the generation of low frequency background firing in a significant subpopulation of DA neurons [7][8][9][10][11][12][13][14][15][16][17][18]51]. However, a number of studies suggest a contribution of Ca 2+ -independent currents to oscillations [26,19,25,27,46]. In accord with these studies, we found that, if considered in isolation, the Ca 2+ -K + oscillatory mechanism provides type II excitability, which is incompatible with the experiments reproduced above. In order to study this, we block the subthreshold sodium current to isolate only the Ca 2+ -and Ca 2+ -dependent K + currents that constitute the oscillatory mechanism. The reduced model is described by the following system of equations: The interaction of the L-type voltage-gated Ca 2+ and the calcium-dependent potassium currents periodically brings the neuron to the spike threshold and generates a metronomic firing activity pattern. Application of an inhibitory input (GABA) to this model reduces the amplitude of voltage oscillations instead of decreasing the firing frequency gradually (Fig 3). This transition is typical for systems where the oscillations are terminated via an Andronov-Hopf bifurcation: the oscillatory trajectory (limit cycle) decreases in amplitude and merges with an equilibrium state. An abrupt transition to zero frequency quiescence upon hyperpolarization suggests a type II excitability of the neuron. Thus, these currents are important for producing metronomic firing, but other intrinsic currents allow for a gradual frequency decrease during application of a hyperpolarizing input observed experimentally.
The role of a subthreshold sodium current. We found that the subthreshold sodium current is necessary for gradual frequency decreases during application of a hyperpolarizing input, as well for a frequency range that spans the observed control frequencies of the DA neurons during a balanced tonic NMDAR and GABAR activation. The reduced firing frequency in the balanced state comes about because the inputs create a slow "bottleneck" effect in the subthreshold voltage range, where the hyperpolarizing inputs nearly cancel the depolarizing ones (Fig 4, also see methods, subthreshold currents section for a more detailed description) and expand the interspike interval. This input balance is achieved due to the contribution of the subthreshold sodium current into the pacemaking mechanism of the DA neuron. By contrast, in the reduced model that includes only Ca 2+ and Ca 2+ -dependent K + currents into the mechanism, the inhibitory input cannot restore appropriate frequency, but instead blocks the voltage oscillations. The inclusion of the subthreshold sodium current allows the firing frequency to vary without compromising the oscillatory mechanism. It allows us to reduce the frequency to an arbitrary low value upon application of a hyperpolarizing input and thus leads to type I excitability of the DA neuron.
A mechanism for NMDA-GABA balance: SNIC bifurcation. Mathematical analysis allows for a better explanation of how the subthreshold sodium current changes the response of the neuron to a combination of excitatory and inhibitory inputs. First, we further reduce the model by separating the slow dynamics from the fast dynamics and removing the fast sodium and the delayed rectifier potassium spike-producing currents. We now define the model as having fired a spike when the voltage crosses a putative spike-threshold (set at -40 mV).  Application of GABA leads to termination of oscillations through Andronov-Hopf bifurcation in a DA neuron, which pacemaking produced by pure Ca 2+ -K + oscillatory mechanism. (A) Two-parameter bifurcation diagram showing transition from type II to type I excitability produced by increasing subthreshold sodium current. (B) F-I curves for two cases shown in (C1) and (C2). Subthreshold sodium current enables smooth frequency decrease to zero upon GABAR activation. One-parameter bifurcation diagrams in cases when subthreshold sodium current is absent (C1) and present (C2). SNIC stands for a saddle node bifurcation on invariant circle, B-T is Bogdanov-Takens bifurcation, A-H is Andronov-Hopf bifurcation, LC stands for a limit cycle. the experiments, in which blockade of the spike-producing currents do not significantly change either the frequency of background oscillations [4] or the NMDA-evoked high frequencies [40]. Additionally, model behavior with fast sodium current and when it is blocked is shown in S2 Fig. It illustrates that full DA neuron is type I excitable as well as the reduced model (without spike-producing currents, but with the subthreshold Na + current). The reduction to the slow-dynamics subsystem decreases the number of variables in the model and enables standard nullcline analysis because it applies only to two-dimensional systems [29,52] (see methods for the description of the nullcline analysis).
Our reduced Ca 2+ -K + model for voltage oscillations (e.g. [4]) shows that the transition to quiescence at a hyperpolarized voltage occurs via an Andronov-Hopf bifurcation (Fig 3), in which oscillations disappear without a decrease in the frequency. Fig 6 presents standard nullcline analysis of the model, which explains oscillatory behavior and bifurcations mechanistically. In Fig 6A, the Andronov-Hopf bifurcation occurs as the voltage nullcline shifts down and simultaneously to the right, so that its intersection with the Ca 2+ nullcline moves across the minimum. In the model with the subthreshold sodium current and Ca 2+ leak current (Fig 6B), the minimum of the voltage nullcline is further away from the steep part of the Ca 2+ nullcline, so that, when the voltage nullcline shifts down, its minimum touches the flat part of the Ca 2+ nullcline. The proximity of the minimum of the voltage nullcline and the bottom part of the Ca 2+ nullcline creates a "bottleneck" effect: The closer the nullclines, the smaller the vector field (the rate of change) in this neighborhood. The limit cycle is channeled through the gap between the nullclines and, accordingly, the oscillation evolves slowly. In the limiting case, when the minimum of the voltage nullcline touches the bottom part of the Ca 2+ nullcline, a saddle-node on invariant circle (SNIC) bifurcation occurs: two equilibrium states, a stable (node) and an unstable (saddle) emerge, interrupt the limit cycle, and the period becomes infinite. The closer the bifurcation parameter g GABA to the bifurcation value, the more time the voltage spends in the bottleneck, creating a long interspike interval. Thus, by introducing the subthreshold sodium current, we change the bifurcation that leads to the quiescence at hyperpolarized potentials from Andronov-Hopf to SNIC.
The role of Ih current. DA neurons are often identified by a prominent hyperpolarization-activated cation current (Ih), which gives a voltage "sag" upon injection of a hyperpolarizing current. It has been shown that Ih-expressing DA neurons exhibit smooth frequency decrease, pointing to a type I excitability [35,36]. However, the excitability type has been shown to vary depending on the intrinsic currents and network connectivity [42,43]. For The subthreshold sodium current changes the type of transition to hyperpolarized rest state induced by GABAR activation. In both cases the growing GABAR conductance leads to a downward shift of the voltage nullcline (solid folded curve). (A) In the pure Ca 2+ -K + mechanism for voltage oscillations, inhibition leads to a transition to the rest state through an Andronov-Hopf bifurcation, which occurs with little change in the firing frequency. The Andronov-Hopf bifurcation is defined as the disappearance of a closed trajectory representing firing (limit cycle) by shrinking in amplitude and merging with an equilibrium state. (B) If the Ca 2+ -K + mechanism is augmented by a subthreshold Na + current, the transition occurs through a Saddle-Node on Invariant Circle bifurcation (SNIC), which corresponds to a gradual decrease in the frequency to zero. The SNIC bifurcation is defined as the emergence of new equilibrium states that interrupt the limit cycle.
doi:10.1371/journal.pcbi.1005233.g006 example, modeling results suggest that changes in the intrinsic currents, e.g. L-type Ca 2+ current, can switch the excitability type of the DA neuron [44,45]. Our model suggests that DA neurons display type II excitability in the presence of the Ih current. However, at small values of Ih conductance (gh = 1-6 mS/cm^2), the model behaves very similarly to the one without the Ih current, which might be taken for type I excitability in experiments. For small values of Ih conductance, the firing frequency of the simulated DA neuron decreases abruptly with the increase in the hyperpolarizing current. However, the discontinuity occurs at very low frequencies (see Fig 7), which might appear still as a gradual frequency decrease in the experiments. As the strength of Ih conductance increases (gh>6 mS/cm^2) the discontinuity in the F-I curve becomes clearer. For high conductances of Ih, low frequency tonic firing cannot be produced, which might affect the basal DA levels. This apparent switch in the excitability type might be important in the altered states of the DA system, when the Ih current is potentiated, for example, by EtOH [53,54].
As we describe earlier, the smooth frequency decrease and the transition to the rest state upon application of hyperpolarizing current or activation of GABAR occurs via a SNIC bifurcation ( Fig 7C3). As soon as we include the Ih current into the model, the SNIC bifurcation splits into a saddle-node bifurcation of limit cycles (LC Fold) and two subcritical Andronov-Hopf bifurcations as shown in a two-parameter bifurcation diagram in fig 7A (see caption for the definitions of bifurcations). Thus, the excitability type changes to type II as we introduce Ih. However, for small values of its conductance gh, the two bifurcations are very close together and nearly indistinguishable, thus, the discontinuity in the F-I curve is in a very narrow range close to zero Hz, producing an illusion of a smooth transition to the rest state.
The bifurcation scenario for the values of gh less than 7mS/cm^2 is complex, although this complexity does not affect experimental observations as the affected limit cycle is unstable. As the magnitude of the negative applied current increases, the unstable limit cycle emerging from the Andronov-Hopf bifurcation, disappears as it collides with a saddle equilibrium state in a homoclinic bifurcation (not marked) and then appears again as the saddle disappears (SN1 in fig 7). The amplitude of the unstable limit cycle grows with the further increase in the negative applied current and, finally, it merges with the stable limit cycle for spiking and annihilates it at the fold bifurcation for limit cycles (LC Fold, Fig 7C2). For higher values of gh>7mS/cm^2 the curve of equilibrium states becomes monotonic, the saddle-node bifurcations of equilibrium states disappear and the transition from spiking to the rest state occurs via the LC fold bifurcation (Fig 7C1). Fig 7 illustrates representative examples of one-parameter bifurcation diagrams and F-I curves for all three cases: gh = 0, 0<gh<7 and gh> = 7.
For a range of negative applied currents, a stable limit cycle coexists with a stable equilibrium state, creating bistability. Thus, depending on the initial conditions, the neuron will be either in a rest (at equilibrium point) or in a repetitive spiking state (at limit cycle). The bistability region grows with the increase in the conductance of the Ih current ( Fig 7A). Variations in current strength back and forth across this range will cause the neuron to jump from one state to the other. To check predictions generated by our model regarding the presence of hysteresis with respect to the applied current, slowly rising and falling current ramps can be applied to the DA neuron to switch it from the rest state to the repetitive firing mode and vice versa. The switch points should be different for falling and rising stimuli as shown in fig 7B. Despite the range of current intensities for coexistence is small relative to the total range for repetitive firing, it affects neuron behavior near the threshold. For example, once the current intensity reaches the value required for the onset of repetitive spiking, small perturbations of current will not switch the neuron back to the rest state. Thus, the presence of hysteresis makes spiking near the threshold more robust.

Changes in the type of excitability caused by synaptic inputs
In vivo, the type of excitability may change due to tonic synaptic inputs [42], and next we explore how this change occurs in the DA neurons. AMPA receptor may co-activate together with NMDA and GABA receptors in vivo. By contrast to NMDAR, conductance of which is voltage-dependent, AMPA and GABA receptor currents are purely ohmic. Their combination LC Fold is a saddle-node bifurcation of limit cycles, in which two oscillatory solutions, stable and unstable, emerge. A-H is the Andronov-Hopf bifurcation, in which a limit cycle shrinks in amplitude and merges with an equilibrium state. SN is a saddle-node bifurcation of equilibrium states, in which two equilibria (rest states) emerge. SNIC is a saddle-node bifurcation of equilibria on invariant circle, which is the same SN bifurcation that occurs on a limit cycle and, therefore, interrupts it. (B) F-I curves for three values of gh: gh = 0mS/cm^2, gh = 3.5mS/cm^2 and gh = 7mS/cm^2. (C1-3) Representative one-parameter bifurcation diagrams of three different bifurcation scenarios occurring at gh = 0mS/cm^2, 0<gh<7mS/cm^2 and gh> = 7mS/cm^2. Black curves represent equilibria. Red curves represent minima and maxima along a limit cycle. Solid stands for stable and dashed for unstable solutions. Dotted arrows represent bifurcation transitions. In C1, spiking is blocked as Iapp decreases through the LC Fold bifurcation. If Iapp increases through the same interval, the LC Fold bifurcation stays unnoticed because the equilibrium state remains stable, and only after the A-H bifurcation, spiking emerges. This is called hysteresis, and it creates bistability, in which either spiking or the rest state can be observed depending on the initial conditions for the voltage and Ca 2+ concentration. In C2, the curve of equilibrium states folds and two saddle-node bifurcations of equilibria occur. The equilibria interrupt the unstable limit cycle (homoclinic bifurcation), but do not affect the stable limit cycle yet. In C3, the equilibria interrupt the stable limit cycle and destroy it in the SNIC bifurcation. is also an ohmic current: Where g eff = g AMPA + g GABA is a combined synaptic conductance, and is a synaptic reversal potential. Fig 8 shows the frequency distribution and the type of bifurcation at the transition to the rest state on the plane of these two parameters: conductance and the reversal potential of the ohmic synaptic current. For instance, if the AMPA receptor is blocked, the reversal potential coincides with the GABAR reversal potential, which is in the range of -90 mV to -70 mV [55]. In this range, an increase in the conductance E eff leads to a gradual decrease in the frequency to zero and a transition to the rest state via a SNIC bifurcation. By contrast, at higher reversal potentials, the frequency drops to zero abruptly and the transition corresponds to an Andronov-Hopf bifurcation. This suggests a transition to type II excitability for the DA neuron. Thus, elevated GABAR reversal potential or tonic activation of AMPAR leads to a switch in the excitability to type II.
Activation of GABAR cannot rescue firing blocked by tonic AMPA AMPAR-mediated input induces depolarization block in the DA neurons and firing cannot be restored. This is consistent with the previous modeling studies showing that NMDA can elicit bursting [56][57][58] or burst envelope [59], while AMPA abolishes high frequency firing. The dynamical explanation is that AMPAR activation shifts the minimum of the voltage nullcline across the Ca 2+ nullcline, so that for high AMPAR conductance values (as well as positive applied currents), voltage oscillations decrease in amplitude and depolarization block occurs  [37,38,40]), Application of GABAR-mediated input shifts the voltage nullcline even further to the right and makes the stable equilibrium more robust. Therefore, the region of parameters for which spiking is produced is much smaller for combined application of AMPA and GABA than for combined NMDA and GABA activation (compare Figs 2B and 9A). Therefore, the prediction of our model is that a disinhibition burst can be supported by tonic background activation of NMDA but not AMPA receptors. Please note that nullclines can be produced only for the model without the fast sodium current (see S3 Fig for  Activation of NMDAR shifts the boundary between the excitability types to lower values of the GABA reversal potential Together with AMPA and GABA receptors, the NMDAR may be also co-activated, since glutamate binds to both AMPA and NMDA receptors. To make the analysis of the excitability type possible in the parameter space of all three synaptic currents, we further perform 2-dimensional bifurcation analysis. The point marked Bogdanov-Takens bifurcation in Fig 8 is a good predictor of the type of excitability at the boundary between spiking and the rest state. Mathematically, it is defined as a junction of the SNIC bifurcation and the Andronov-Hopf bifurcation, as it appears in the figure. In Fig 10, we plot this point as a function of the NMDA  Fig 10 shows that the separation between the types of excitability shifts to lower values of the combined reversal potential for the AMPAR and GABAR currents as the NMDAR conductance increases. However, this shift quickly saturates and is restricted to the range of GABAR reversal potentials. Thus, similar to Ih, NMDAR activation may switch the type of excitability of the DA neuron from type I to type II in a certain window of other synaptic currents received by the neuron.

DA release and synchronization in a population of heterogeneous DA neurons is influenced by tonic background synaptic currents
To illustrate the importance of changing DA neuron type of excitability, we simulated heterogeneous populations of DA neurons under two conditions: 1) in control (in the absence of the tonic synaptic inputs), and 2) during the tonic influence of AMPA and GABA inputs. The DA population is electrophysiologically heterogeneous [60][61][62], and its uncoordinated activity produces a homogeneous low-level DA concentration. In order to have similar firing rates and basal DA levels in both cases, we balanced the increase in firing rate produced by the application of AMPAR conductance with GABAR conductance (note that GABA can balance AMPA only for a very limited range of values). In both cases, DA neurons received correlated fluctuating NMDA inputs (Fig 11A, see methods for the detailed description of the inputs). Our simulations show that the population of DA neurons that receive the background synaptic tone produces higher DA levels in response to bursty correlated NMDA The boundary between type I and type II excitability in the DA neuron shifts down with respect to the reversal potential of an ohmic synaptic current as NMDAR conductance grows. (B1-2) Two voltage traces illustrating the excitability types. For any synaptic reversal potential E eff and NMDAR conductance from the grey region of type I excitability on panel (A), ramping the ohmic synaptic conductance g eff increases the interspike intervals and then blocks firing. For the parameters from type II region, ramping the ohmic synaptic conductance blocks firing without the decrease in the frequency.  input than a population without the synaptic tone ( Fig 11D). As described above, DA neurons are type I excitable in the absence of synaptic tone, while AMPAR activation switches DA neuron excitability to type II. Thus, the transition from type I to type II excitability of the DA neurons is accompanied by higher dopamine release in response to a correlated synaptic input. The higher responsiveness is partially due to a greater synchronization of the DA neurons receiving the synaptic tone, as evident by the higher number of peaks in their summed activity in Fig 11C. Type II neurons display more robust synchrony when they receive a common input, even in the presence of independent noise [33,63]. Thus, in vivo background synaptic tone might be important not only for regulating basal DA levels, but also for the responsiveness of the DA neurons, so that they are more ready to produce coincident bursts in response to correlated synaptic inputs.

Discussion
The excitability type of the DA neuron The type of excitability for a neuron determines the neurons' responses to stimuli and their dynamics in a population (phase response curve, synchronization, resonators vs. integrators). The type of bifurcation determines the type of excitability: neural oscillations that arise via an Andronov-Hopf bifurcation have type II excitability, while those appearing via SNIC have type I excitability [52,65]. Based on the bifurcation analysis and frequency responses to hyperpolarizing inputs (negative injected current and GABAR conductance), we have shown that in control conditions, the simulated DA neuron is type I-excitable. It's known that the type of excitability in vivo may be different from in vitro [42]. In vivo, DA neurons display irregular low-frequency activity occasionally interrupted by high-frequency bursts. This low-frequency regime may reflect the balance of inhibitory and excitatory inputs. We found that, in the most prominent low-frequency regime, the DA neuron is type I excitable, in either high or low synaptic conductance states.
The baseline level of dopamine is important for the normal function of the brain. The level is determined by the background activity of the DA neurons. This activity is intrinsically generated by the neuron and controlled by its synaptic inputs [66], reviewed in [67]. Thus, the capacity of the DA neuron to adjust its firing rate according to the inputs is vital. The graded response curve of a type I-excitable neuron, as opposed to an on-off response of a type II neuron, provides this capacity. Accordingly, at every level of excitation provided by NMDAR input, inhibitory GABA synaptic conductance can balance it and bring the frequency down to an arbitrary low value. A similar hyperpolarizing current activated by dopamine D2 receptors on the DA neuron has also been shown to slow down its firing rather than abruptly block it all together [68]. These are very important autoregulatory functions of the DA system that allow it to adjust basal DA levels in target areas.

The role of the subthreshold sodium current
Persistent sodium current is known to amplify subthreshold oscillations [69] and increases neural excitability of DA neurons by contributing to spontaneous depolarization in between the spikes [25]. Consistent with experimental observations, the subthreshold sodium current population receiving the tonic synaptic inputs. This is also illustrated in the insert zoomed in on the neural activity at the beginning of the burst. (D) Cumulative synaptic DA concentration produced by activity of the DA population with (red) and without (black) tonic background synaptic inputs. The DA neuron population that recieves tonic inputs releases more dopamine in response to NMDAR stimulation. increases the firing rate of the DA neuron in the model. Additionally, we found that the current is necessary for achieving gradual frequency decrease upon application of hyperpolarizing current, thus, maintaining type I excitability of the DA neuron. The type of excitability is determined by the internal properties of the currents contributing to pacemaking in the DA neuron. L-type Ca 2+ and SK-type Ca 2+ -dependent K + currents are the core currents that traditionally constitute this mechanism [7][8][9][10][11][12][13][14][15][16][17][18] (but see the section on the mechanisms of DA neuron pacemaking). However, our model shows that the mechanism results in type II excitability, in which a hyperpolarizing current blocks the voltage oscillations without restoring a low frequency. Our explanation is that, without the contribution of further currents, the steep part of the Ca 2+ nullcline is very close to the minimum of the voltage nullcline (Fig 4A), because they reflect the same event: opening of the Ca 2+ channel. This positions the system close to the Andronov-Hopf bifurcation, which occurs whenever the minimum of the voltage nullcline moves across the Ca 2+ nullcline. When the subthreshold sodium current is included into the mechanism, the minimum of the voltage nullcline reflects opening of this current, and it is shifted away from the steep part of the Ca 2+ nullcline (Fig 6B). This shifts the system away from the Andronov-Hopf bifurcation. Now, a downward shift of the voltage nullcline following inhibitory inputs moves its minimum across the flat part of the Ca 2+ nullcline and produces a SNIC bifurcation. In this transition, the firing frequency gradually reduces to zero, and this allows a balance between NMDAR and GABAR conductances and restores the background firing frequency.

Related studies of NMDA-GABA balance in DA neurons
The influence of NDMA and GABA receptor conductances on the DA neuron have been studied in several papers [23,22,70]. The compensatory influence of NDMA and GABA receptor activation on the firing frequency has been predicted in modeling studies by Komendantov et al. [23]. Lobb et al. [70] modified our previous model [24] to capture the balance of the inhibition and excitation and disinhibition bursts. In these models, a high maximal frequency (>10 Hz) can be achieved by tonic activation of the AMPAR or in response to depolarizing current injection, which contradicts experimental observations [7,34]. To the contrary, a number of experimental studies suggest that stimulation of NMDA receptors evokes a burst of high-frequency firing, whereas AMPA receptor activation evokes modest increases in firing [37][38][39][40] (but see [61,71]). This is an important distinction, which impacts the excitability of the neuron. Further, in Komendantov et al. [23] and Canavier & Landry [22], the NMDAR conductances were restricted to dendrites, whereas GABAR conductance was somatic. The mechanism of frequency rise during dendritic application of NMDA is different from the mechanism of response to somatic NMDAR stimulation [56,24]. Somatic NMDAR stimulation has been shown to elicit high-frequency firing in earlier experiments [46,40] and used to achieve the NMDA-GABA balance [35]. Here, we base a new model on our previous model [56] that presented a mechanism for somatically-induced high-frequency firing in a reconstructed morphology first and reduced it to a single compartment. In the current model, we have integrated the mechanism for high-frequency firing together with the balance of NMDAR and GABAR activation.

Mechanisms of DA neuron pacemaking
The mechanism of low frequency pacemaking in the DA neurons has been extensively studied. However, it is still a matter of on-going debate in the literature since different experimental results lead to contradicting conclusions, proposing that different currents are critical for the DA neuron spontaneous firing. In a number of experimental and modeling studies it has been shown that spontaneous tonic firing relies on the interactions between the voltage gated calcium (Ca 2+ ) and SK-type Ca 2+ -dependent potassium (K + ) currents [7][8][9][10][11][12][13][14][15][16][17][18][19]. Wilson and Callaway [4] and later Chan et al. [19] showed that calcium-driven slow oscillatory potentials (SOPs) drive the spiking rate of the SNc DA neurons. Chan et al 2007 [19] also showed that dependence of pacemaking on Ca 2+ oscillations changes with the age. Particularly, TTX blocks slow oscillations in juvenile neurons, but not in adult neurons, which is related to the change in density of Ca 2+ channels with age. Ping and Shepard 1996 showed that the frequency of SOPs after the application of TTX is approximately the same as the frequency of spiking. In contrast, in a more recent study, Guzman et al. [26] demonstrated that, in a number of DA neurons, SOPs and spiking frequencies are weakly correlated, and that TTX inhibits spontaneous oscillatory potentials pointing to the importance of sodium currents for pacemaking. A number of other studies also suggest that sodium channels are highly involved in controlling spontaneous DA neuron frequency [26,25], especially in the VTA DA neurons [27].
The sources of the apparent discrepancy in the experimental results were investigated by Drion and colleagues [21]. Based on the combination of experimental and modeling approaches, authors suggested that calcium and sodium currents likely cooperate to produce pacemaking and prevailing mechanism depends on the density of the ion channels in the neuron. Further, authors showed that the lack of correlation between spikes and SOPs does not lead to a conclusion that generating mechanisms are different. The complementary role of the two mechanisms, notably if they co-exist in the same cell or represent pacemaking in distinct populations, is a matter of on-going research in the field.
For the sodium-based pacemaking mechanism, it is still unclear what hyperpolarizing current provides the long interspike interval when the SK current is not functional. This renders the Ca 2+ -independent oscillatory mechanism incomplete and is a matter of a future investigation. Control of the firing by the SK current has been shown to be stronger in the DA neurons positioned more laterally in the midbrain [15]. Thus, we focus on the subpopulation of DA neurons that are more abundant in the substantia nigra pars compacta (SNc) than the ventral tegmental area (VTA).

Intrinsic Ih current and synaptic inputs can switch the excitability type of the DA neuron
Changes in intrinsic currents can affect the excitability type and, thus, computational properties of the DA neuron. For instance, we observed that potentiation of Ih current promotes type II excitability of the simulated DA neuron (Fig 7). Further, we show that Ih current can induce bistability in the DA neuron, and the bistability region increases with the increase in Ih conductance (Fig 7B). This can affect the behavior of the neuron near the boundary between spiking/resting states, particularly, if the current reached the value necessary to induce spiking onset, a small perturbation in the current will not silence the neuron. This makes spiking more robust near the threshold. In addition to the contribution of Ih current to pacemaker activity, has been shown in DA neurons [72], as well as in the other neuronal types that Ih induces intrinsic subthreshold resonance [73][74][75]. Thus, augmentation of Ih current increases oscillatory behavior of the DA neurons, as well as their synchronization in response to excitatory pulses. However, low-frequency tonic firing could not be maintained at high conductances of Ih current, likely affecting background DA levels.
Further, the influence of tonic synaptic inputs can also change the transition to the rest state and, therefore, be described as altered excitability. Tonic activation of AMPA receptors or an elevated reversal potential of the GABAR conductance may make the low-frequency balanced state unreachable. The reason for that is a transition to type II excitability: firing is blocked at higher frequencies. A model prediction that follows from this result is that tonic AMPAR activation induces depolarization block and firing cannot be rescued by application of GABA. Our explanation is that shunting is so strong that opening of the subthreshold sodium current cannot sustain the voltage growth. In other words, these changes unfold the voltage nullcline and bring its minimum close to the steep part of the Ca 2+ nullcline (Fig 9B). This primes the system for the Andronov-Hopf bifurcation responsible for type II excitability. Further, we found that NMDAR activation also biases the neuron towards type II excitability (Fig 10). Although the type may change as parameters shift away from the boundaries of the firing region [76], together, these results suggest that in high-frequency regimes the DA neuron displays type II excitability. This switch in excitability type may play a role in a transient increase in DA concentration in response to salient stimuli as it is easier to synchronize type II neurons by an excitatory input. Fig 11 supports this hypothesis by showing higher transient DA release produced by heterogeneous population of DA neurons receiving synaptic AMPA and GABA tones than in the absence of synaptic tone. Thus, correlated excitatory synaptic inputs are more likely to evoke robust coincidence DA release when DA neurons display type II excitability.

Implications of the DA neuron excitability type on computations and possibly different behavioral states
A growing body of literature links the type of excitability to neural coding [29,31,41,65,[77][78][79]. For instance, Prescott et al. [80] suggested that type I neurons are best suited for coding stimulus intensity. Hence, the DA neuron is designed for encoding the intensity of the tonic depolarizing and hyperpolarizing inputs by its smooth frequency dependence. This further supports and augments a recently found unique computational role for the DA neuron: it performs subtraction of inhibitory and excitatory inputs [81]. The operation is optimal to calculate unpredicted value of an event but rarely observed and hard to implement in the brain. The first type of DA neuron excitability is necessary to quantitatively encode the level of input by the firing rate and perform the subtraction. Several studies, for example Eshel et al. 2016 [82] and Tobler et al. 2005 [83], showed that activation of DA neurons gradually increases with the increase in the reward value. We attempted to reproduce this experimental result in our model. For the simulations, we assumed that the reward value is proportional to the strength of NMDA input coming to the DA neurons (g NMDA ). GABA inputs do not seem to be a plausible candidate because Eshel et al. 2016 [81] showed that firing of the GABA neurons, as opposed to the DA neurons, does not vary consistently with the reward value. We calculated firing frequency dependence of the type I/ type II DA neurons on the input strength (Fig 12). Type II DA neurons were modeled by applying tonic AMPA along with NMDA conductance (additional excitation was compensated by increasing GABA conductance). Type II DA neurons are unable to encode low reward values as their firing abruptly drops to zero with the decrease in the input strength. In contrast, frequency dependence of type I DA neurons resembles the experimental curve of the DA neuron frequency dependence on the reward value shown in [82]. The difference can also be seen in the raster (Fig 12A and 12C) and frequency responses (Fig 12B and 12D) to a transient increase in the input strength, representing the reward value. First peak in the frequency response represents salience, and was simulated by applying constant level of NMDA for all the values, while second peak represents the reward value, and was simulated by scaling NMDA input accordingly. There is a gap in a frequency response to low reward values of type II DA neurons (Fig 12D). Thus, type I DA neurons best encode the value of an event, i.e. the difference between predicted and received reward.
On the other hand, spike timing of type I neurons in response to weak or noisy transient inputs is not reliable [84,85], while temporal precision of type II neurons is much higher. The ability of the DA neuron to switch excitability type from type I to type II under certain synaptic inputs might play a significant role in producing enhanced transient DA release, since it likely relies on the precise coordinated activity of the DA neurons. Multiple drugs of abuse, including EtOH (e.g. [86]) evoke transient increases in the DA concentration in nucleus accumbens.
Using our model, we show that EtOH shifts excitability of the DA neurons to type II and induces higher DA release. EtOH acts on multiple intrinsic and synaptic currents. Mainly, it enhances I h current [87,88] and increases AMPA/NMDA ratio [89]. Moreover, it increases GABA release onto DA neurons [90]. Thus, we modeled EtOH action by increasing I h , AMPA and GABA receptor currents. Ih and AMPAR currents switch DA neuron excitability to type II and, therefore, promote synchronization in the population of DA neurons in response to noisy excitatory inputs (Fig 13). This is one of the mechanisms that can produce higher DA transients. By switching DA neuron to type II, EtOH can increase the motivational potential of a stimulus because the same excitatory input produces enhanced DA signal under EtOH. In other words, neutral stimulus can become salient after EtOH exposure. Our modeling prediction regarding the excitability switch after EtOH could be tested in-vitro by studying how the shape of F-I curve changes after EtOH. In in-vivo conditions it could be checked by stimulating DA neurons before and after EtOH exposure with a chirp pattern signal in order to check for spiking resonance. To test whether DA neurons better synchronize after EtOH, a phase response curve (PRC) or a spike triggered average (STA) could be calculated. Presence of the negative PRC component and narrow STA indicates that neurons are more amenable to synchronization by common synaptic noise.
In conclusion, our results predict that DA neurons can exhibit traits of both integrators and resonators and these traits are modulated by intrinsic and synaptic conductances. Depending on the current constitution, DA neurons can perform rate coding by integrating slow variations in the inputs and adjust basal DA concentration or they can detect transient coherent changes in the inputs and synchronize for producing robust DA transients.

Methods
The biophysical model of the DA neuron is a conductance-based one-compartmental model modified from [91] where v is the voltage and c m is the membrane capacitance. There are eight intrinsic currents of the DA neuron: a calcium current (I Ca ), a calcium-dependent potassium current (I KCa ), a potassium current (I K ), a direct rectifier current (I DR ), a subthreshold sodium current (I sNa ), a hyperpolarization-activated current (I h ), a fast sodium current (I Na ), and a leak current (I leak ). The first subgroup of intrinsic currents: I Ca , I KCa , I K , IsNa and I h constitute a pacemaking mechanism of the DA neuron. The second subgroup of the intrinsic currents (I Na , I DR ) is responsible for spike generation. The last subgroup includes synaptic currents: the excitatory α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic (AMPA) and N-Methyl-D-aspartate (NMDA) receptor currents (I AMPA and, I NMDA respectively) and the inhibitory γ-Aminobutyric acid (GABA) receptor current (I GABA ). Synaptic inputs can produce bursts and pauses in firing.

Intrinsic oscillator
The main currents of the model that produce pacemaking activity of DA neuron are an L-type voltage-dependent calcium current (I Ca ) and an SK-type calcium-dependent potassium current (I KCa ). Gating of the calcium current is instantaneous and described by the function: Calibration of the calcium gating function reflects an activation threshold of an L-type current, which is significantly lower in DA cells than in other neurons (~-50mV; [4]). Calcium enters the cell predominantly via the L-type calcium channel. Contribution due to the NMDA channel is minor [92]. Thus, calcium concentration varies according to the second equation of the system (1). It represents balance between Ca 2+ entry via the L channel and a Ca 2+ component of the leak current, and Ca 2+ removal via a pump. In the calcium equation, β is the calcium buffering coefficient, i.e. the ratio of free to total calcium, r is the radius of the compartment, z is the valence of calcium, and F is Faraday's constant. P ca represents the maximum rate of calcium removal through the pump. A large influx of Ca 2+ leads to activation of the SK current, which contributes to repolarization as well as afterhyperpolarization of the DA cell. Dependence of the SK current (I KCa ) on calcium concentration is modeled as follows: The neuron is repolarized by the activation of a large family of voltage-gated potassium channels. In addition to the already described potassium current, the model contains voltagedependent K current (I K ). Conductance of this current is given by a Boltzmann function: The DA neuron expresses voltage-gated sodium channels that carry a large transient current during action potentials (a spike-producing sodium current) and a non-inactivating current present at subthreshold voltages (a subthreshold sodium current). Even though the persistent subthreshold sodium current is much smaller than the transient spike-producing current, it influences the firing pattern and the frequency of the DA neuron by contributing to depolarization below the spike threshold [25]. We modeled the voltage dependence of the subthreshold sodium current as follows: The kinetics and the voltage dependence of the subthreshold sodium current were taken from [93].
The majority of DA neurons express a hyperpolarization-activated nucleotide-gated (HCN) inward cation current (I h ). The HCN current contributes to spontaneous firing of subpopulations of DA neurons [94]. The activation variable of I h is governed by a first-order ordinary differential equation (the third equation of the system (1)).
The maximal activation of I h current is described by the following voltage-dependent equation [88] The voltage-dependent time constant is described by The leak current (I leak ) in the model has the reversal potential of -35 mV, which is higher than in the majority of neuron types. In DA neurons, several types of depolarizing, nonselective cation currents are expressed, which likely contribute to depolarization during interspike intervals.

Model calibration
Subthreshold currents. Calibration of the subthreshold Ca 2+ -K + oscillatory mechanism remained very similar to previously described [89]. Experiments show that subthreshold sodium currents contribute to depolarization towards the spike threshold [25,27]. Accordingly, the addition of the subthreshold sodium current into the model increases the background firing frequency of the DA neuron. However, in the model this leads to a significant increase in maximal frequencies achieved during application of a constant depolarizing current or AMPA, whereas in experiments these frequencies have been shown to be limited by approximately 10 Hz [34]. We preserve this limit in the model by accounting for a Ca 2+ component of the leak current. Addition of the small Ca 2+ leak current moves the flat part of the Ca 2+ nullcline up and produces a slight upward bent of the Ca 2+ nullcline, which was not essential for the SNIC bifurcation. This reinforced the negative feedback loop through the Ca 2+ -dependent K + current and limited the maximal frequencies. NMDAR activation increases the frequency, and our calibration of the model allowed an inhibitory or hyperpolarizing input to reduce the frequency to an arbitrary low value. Ca 2+ entry through NMDA receptor was omitted to implement the spatial segregation of NMDAR and SK channels as in the previous models [56,23].
Spike-producing currents. We included spike-producing sodium current with maximal conductance g Na ðvÞ ¼ g Na m 3 h and potassium delayed rectifier current with maximal conductance g DR ¼ g DR n 4 . The activation of the Na + current is assumed to be instantaneous and is described by the function The inactivation kinetics is described by the equation for h in System (1), where The delayed rectifier activation variable n is described by the 5th equation in System (1), where a n ¼ À 0:0032ðv þ 5Þ The currents are calibrated to produce a spike per each maximum of the voltage oscillations produced by the pacemaking mechanism without significantly changing the firing rate or pattern.

Synaptic inputs
DA neurons receive glutamatergic (Glu) excitatory drive through AMPA and NMDA receptors and inhibitory drive through GABA receptors. Changes in the membrane potential induced by synaptic conductances are described by the following equation where g NMDA (v), g AMPA , g GABA are the maximal conductances of NMDA, AMPA and GABA receptor currents accordingly, s NMDA , s AMPA , s GABA are gating variables that depend on the input spike trains. The AMPA and GABA conductances are voltage-independent, but the NMDA conductance has voltage sensitivity as in [95] where [Mg 2+ ] denotes the amount of magnesium, taken to be 0.5μM. The low slope of the voltage dependence (m e = 0.062) is critical for the increase in the frequency of spikes or subthreshold oscillations during NMDA application [56].
Modeling of Glu and GABA asynchronous inputs to the DA neuron Glu input. We started our study with investigating the influence of asynchronous synaptic inputs. Asynchronous Glu input to the DA neuron was produced by 35 Poisson distributed spike trains with frequencies of approximately 10 Hz. The number of spike trains was chosen to produce a relatively constant level of NMDA receptor activation, and, at the same time, take into account the effects of convergent synaptic inputs on the DA neuron, by thresholding NMDAR to activate only by coincidence of two or more spikes. The activation of the receptors in response to a synaptic input is described by the following equation where j denotes a dimensionless synaptic input. It is normalized to change from 0 to 1 for 1 ms interval to mimic a single spike in the input. i denotes a receptor type, AMPA or NMDA. GABA input. Asynchronous inhibitory input was produced by a population of 30 GABA neurons. The number of GABA neurons projecting to a DA neuron is not known. To choose this number, we note first that the percentage of GABAergic neurons varies between 12 and 45% in different sub regions of the ventral tegmental area (VTA) [96]. Since the GABA neurons powerfully modulate DA neuron activity via direct, monosynaptic inhibitory connections [55,[97][98][99][100], one can expect multiple GABA neurons to make connections with a single DA neuron. Second, we found that our results are valid for a wide range of the number of GABA neurons and start to change only if the number becomes small (10 or lower) [91].
Voltage dynamics of each GABA neuron is described by the Wang-Buszaki equations of a fast spiking neuron [101], where v i is a voltage of i th GABA neuron in a population. The activation variable m is assumed fast and substituted by its steady-state function m 1  . The activation variable n obeys the following equation, dn dt ¼ a n ðvÞð1 À nÞ À b n ðvÞn; where a n ¼ 0:01ðv i þ29Þ 1À expð À ðv i þ29Þ 10 Þ and b n ¼ 0:0875expð À ðv i þ39Þ

80
Þ. The parameters of these equations were calibrated according to experimental observations [34]. Intrinsic firing frequencies of these neurons were set to a range of 12-22 Hz, similar to what is observed experimentally [34,102,103]. The activity of each GABA neuron contributes to the activation of the GABAR current entering the DA neuron by increasing the gating variable according to the equation where g spike ¼ 1 1þexpð The total receptor activation is a sum of contributions produced by all GABA neurons in a population. A parameter that scales the GABA current is GABAR conductance g GABA . We normalize the GABA gating variable by the number of neurons in order to keep its value in a range between zero and one. Thus, for an asynchronous GABA population, the GABAR will be partially activated and the gating variable will have a low value. Model parameters are given in Table 1.
Modeling tonic synaptic currents. The transition from asynchronous synaptic inputs to tonic synaptic currents was done by reducing NMDAR, AMPAR and GABAR currents to the form of I NMDA = g NMDA (E NMDA − v), I AMPA = g AMPA (E AMPA − v) and I GABA ¼ g GABA ðE GABA À vÞ respectively. This corresponds to setting the receptor activation variables described by eq (13) to 1 to match dynamic clamp experiments [36,50]. Tonic receptor activation (s = 1) is a very good approximation of the asynchronous input ( Compare Figs 1 and 2) Nullcline analysis A multidimensional system can be analyzed with two-dimensional nullcline methods only after its reduction to a two-dimensional system. The description of this method is provided by Strogatz [104] and Rinzel and Ermentrout [29]. The reduction was done by eliminating the oscillations crossed the threshold of -40 mV, as experimentally it was shown that a DA neuron action potential is triggered when the voltage is depolarized to approximately -40 mV [3]. Voltage oscillations that were below these thresholds in the models with and without the spike-producing currents respectively were not counted as spikes and did not contribute to the firing frequency. To analyze firing pattern of simulated DA neuron in the presence of different synaptic currents, we quantified its firing rate and bursting. Mean firing rate of the simulated DA neuron was calculated as an inverse of the mean interspike interval (ISI). To calculate bursting we used ISI coefficient of variation (CV), calculated as the SD/mean of 200 ISIs.