Mechanisms for multiple activity modes of VTA dopamine neurons

Midbrain ventral segmental area (VTA) dopaminergic neurons send numerous projections to cortical and sub-cortical areas, and diffusely release dopamine (DA) to their targets. DA neurons display a range of activity modes that vary in frequency and degree of burst firing. Importantly, DA neuronal bursting is associated with a significantly greater degree of DA release than an equivalent tonic activity pattern. Here, we introduce a single compartmental, conductance-based computational model for DA cell activity that captures the behavior of DA neuronal dynamics and examine the multiple factors that underlie DA firing modes: the strength of the SK conductance, the amount of drive, and GABA inhibition. Our results suggest that neurons with low SK conductance fire in a fast firing mode, are correlated with burst firing, and require higher levels of applied current before undergoing depolarization block. We go on to consider the role of GABAergic inhibition on an ensemble of dynamical classes of DA neurons and find that strong GABA inhibition suppresses burst firing. Our studies suggest differences in the distribution of the SK conductance and GABA inhibition levels may indicate subclasses of DA neurons within the VTA. We further identify, that by considering alternate potassium dynamics, the dynamics display burst patterns that terminate via depolarization block, akin to those observed in vivo in VTA DA neurons and in substantia nigra pars compacta (SNc) DA cell preparations under apamin application. In addition, we consider the generation of transient burst firing events that are NMDA-initiated or elicited by a sudden decrease of GABA inhibition, that is, disinhibition.


Introduction
Dopamine is a key neurotransmitter involved in motivation, memory, and motor function. Two areas in the midbrain are particularly associated with large populations of dopaminergic neurons whose activity transmits dopamine signals: the ventral tegmental area (VTA) and the substantia nigra pars compacta (SNc). The dopamine signal originating from the SNc DA neurons is a necessary component for motor initiation (Burns et al., 1983;Carli et al., 1985;Schultz et al., 1989;Schultz, 2007), and the VTA DA signal is conjectured to convey reward-related information (Corrigall et al., 1992;Schultz, 2002).
The DAergic neurons display a wealth of dynamic behavior: pacemaker-like firing (seen in slice preparations), periods of burst activity, quiescence (due to either strong hyperpolarization or depolarization), and also sporadic activity (see e.g., Grace and Bunney, 1984a,b;Overton and Clark, 1997). In this paper, we examine the mechanisms underlying the range of activity modes found in DA neurons. Given the rich behavior of DA neuronal dynamics and due to their biological importance, a number of modeling efforts have been made to understand DA neural dynamics. Existing models can capture a variety of these behaviors (Wilson and Callaway, 2000;Canavier and Landry, 2006;Kuznetsov et al., 2006), yet are not without their shortcomings. The models of Wilson and Callaway (2000), Canavier and Landry (2006), and Kuznetsov et al. (2006) rely on somatodendritic interactions as underlying burst generation, however, two laboratories have experimental findings that suggest burst firing is somatically induced (Blythe et al., 2009;Deister et al., 2009). These somatically induced bursts have characteristics consistent with in vivo-like bursting, suggesting that a single-compartmental model should be sufficient for generating the observed DA neuronal dynamics.
Here we present a single compartment model of a VTA DA neuron, which captures the essential qualitative behavior of DAergic neuronal dynamics. This revised model improves our earlier work (Oster and Gutkin, 2011) and better captures DA neuronal dynamics. We perform extensive numerical experiments in order to identify conditions for a variety of firing modes and have uncovered novel firing regimes. Other single compartment models have been recently introduced (Drion et al., 2011;Oster and Gutkin, 2011;Ha and Kuznetsov, 2013), and these models differ from the presented model-though do share some common characteristics. For instance, the model of Drion et al. (2011) does not produce endogenous bursting rather burst events are noise driven. Our model dynamics also differ from those of the the dynamics of the Ha and Kuznetsov (2013) model, which predicts a sustained fast firing regime under NMDA-bath-like conditions, which is not inline with the endogenous bursting patterns observed under NMDA-bath application (Johnson et al., 1992). In this work, we consider a wider range of conditions than previous models and uncover dynamics that were heretofore undiscovered by modeling efforts. In this endeavor, we systematically vary the SK conductance and GABA tone to test its effect on firing properties. We also considered alternate sodium activation dynamics and the effects due to bulk NMDA increases and to disinhibition events.
The generation of DA neuronal activity patterns is dependent on a wide number of factors. At the most basic level, the DA pattern is dependent on both the intrinsic properties of the cells and the network activity that innervates the DA cells. Mameli-Engvall et al. recorded and observed an ensemble of activity modes that mouse VTA DAergic neurons express in vivo (Mameli-Engvall et al., 2006), shown in Figures 1A,B. Mameli-Engvall et al. classified neurons by their firing patters as either low or high firing (displaying a frequency of less than or greater than 5 Hz, respectively) and further categorized neurons as either low or high bursting modes (less than 20% being low burst firing and high bursting being greater than 20% of spikes within burst events). Differences in the activity modes are evident when plotting the distributions of the associated interspike intervals, shown in Figure 1B. In addition, Marinelli et al. (2006) present A C B FIGURE 1 | A plot of experimentally observed DA cell firing modes from in vivo rodent VTA adapted from Mameli-Engvall et al. (2006) contrasted with data from computational simulations. In (A), the distribution of in vivo firing modes. Along the horizontal and vertical axes are mean firing frequencies and percentage of spikes within bursts, above 20% was considered a high burst neuron. Neurons firing below 5 Hz were considered to have a low firing rate, and high firing rate was taken to be neurons firing above 5 Hz. In (B), the associated ISI histogram for the data in (A). In (C), we plot the activity modes obtained from considering a family of neurons with varying SK conductances and receiving different levels of drive, I 0 (here GABA is taken to be 0.15). We arrange the data in a similar fashion to Mameli-Engvall et al. (2006), but with the shape of the data points indicating the relative strength of the SK conductance (neurons with weak, moderate and strong SK conductance are denoted via the use of squares, disks, and asterisks, respectively). The entire high firing regime is composed of neurons with low SK conductances and a significant portion of the low firing high burst regime are of low and medium SK conductance neurons. similar histograms for examining VTA DA neural ensembles. In this work, we show that the SK conductance levels act as one of the major determinants of endogenous DA firing modes, Figure 1C.
Dopamine neurons receive glutamatergic input from the prefrontal cortex and a combination of glutamatergic and cholinergic drive from the pontine nuclei. Both VTA and SNc are subject to substantial GABAergic inhibition. In particular, the VTA contains populations of GABAergic neurons that provide local GABAergic inhibition to the VTA DAergic neurons (in addition to other GABAergic inputs from the striatum). The VTA GABAergic populations are subject to similar glutamatergic and cholinergic currents as the DAergic populations, and there is an interplay between the DAergic and GABAergic activities. This balance of excitation to inhibition plays a role in the VTA DA neural response to nicotine stimulation (Mansvelder et al., 2002;Mameli-Engvall et al., 2006;Tolu et al., 2013).
Here in our computational study, we find that the degree of GABAergic tone, together with the glutamatergic activation through the NMDA receptors and the intrinsic cellular properties, plays a strong role in shaping a neuron's firing characteristics. This is not to say that GABA is the most important signal, it appears as though NMDA is likely the principal pathway to initiate burst firing. However, the make-up of the ionic conductances and the degree of the GABAergic tone influence how readily primed a neuron is to burst. So although NMDA can reliably initiate burst firing in DA neurons, the duration of initial burst firing event and the resultant continued firing pattern is dependent on the strength of the DA neuron's SK conductance.

Conductance-based DA Neuronal Model
We present a single compartmental conductance-based model for the dynamics of a DA neuron (a schematic shown in Figure 2) that combines modified conductance mechanisms previously introduced in Amini et al. (1999), Kuznetsov et al. (2006), Canavier et al. (2007, and Oster and Gutkin (2011). Such a model serves to consolidate these approaches to a simplified framework and allows the use of computer simulation to uncover the mechanisms behind the generation of DA activity modes. Parameter values are listed within the text with any unspecified parameters appearing in the Appendix (Table A1).

The Membrane Potential
The dynamics of the membrane potential are described by where I 0 is a generic term meant to incorporate the effects due to either an applied current in vitro or a general tonal dendritic input in vivo. The L-type calcium current, I L Ca , is counterbalanced with the apamin sensitive, calcium dependent SK-type potassium current, I SK K . The SK current is dependent on both the voltage and calcium level, with calcium concentration denoted by u, and the SK current can be parametrically controlled via the parameter χ APA . Two tetrodotoxin (TTX) dependent sodium currents are included that can be parametrically regulated via the parameter χ TTX : (i) a spike-generating sodium current, I Na , modeled in a Morris-Lecar-like fashion with fast activation, m ∞ , and slow inactivation, h, (ii) and a persistent sodium current, I pers Na , suggested to be of particular import to VTA DAergic neurons (Khaliq and Bean, 2010). The DA neurons receive GABAergic inhibition and receive excitatory glutamatergic drive through their AMPA and NMDA receptors. We include both a general and a delayed rectifier potassium currents and a general leak current.
A number of observed currents were not included for simplicity. For instance, previous studies have shown that the hyperpolarization-activated cation current, I h , has little or no effect on VTA pacemaking (Neuhoff et al., 2002;Khaliq and Bean, 2010), and hence we omitted it from this study. As a consequence, FIGURE 2 | Schematic diagram of the single compartmental conductance-based model of DA neuronal activity. The model includes a generalized intrinsic current whose composition could include in vivo cholinergic drive or in vitro applied currents. The model includes ionic currents for sodium, calcium, and potassium. Both spike generating and persistent sodium currents are expressed, which are both selectively blocked by application of TTX. The potassium currents are composed of a general current, a delayed-rectifier current, an apamin-blocked, calcium-dependent SK-type current. Calcium is predominantly trafficked into the intracellular region via the L-type Ca 2+ channel and pumped out of the cell via a pump, with the intracellular calcium undergoing fast buffering. The intracellular calcium modulates the SK potassium current and the interaction between the dominant L-type calcium current and the SK potassium current forms the mechanism for an underlying oscillation. In addition to these currents, DA cells receive excitatory synaptic drive through AMPA and NMDA receptors. In addition, DA cells receive inhibitory input from GABAergic neuronal populations.
our model does not reproduce the characteristic depolarization sag after a hyperpolarization event.

Spike Generating and Voltage Dependent Currents
The fast processes associated with spike generation involve sodium. The sodium dynamics are taken to be Morris-Lecar-like with fast activation and a slow inactivation gating variable h; so that with the conductanceḡ Na = 150 and the reversal potential of E Na = 55 mV. The fast activation is taken to be of the form with p 2 = −14 mV, p 3 = 11.9, and the inactivation gating variable evolves according to In addition, we include a delayed-rectifier potassium current of the form where the conductance is given byḡ K = 4 and the reversal potential E K = −90 mV. The associated gating variable evolves according to with α n (V) = n a1 2 1 − tanh n a2 − V n a3 and β n (V) = n b1 2 1 + tanh n b2 − V n b3 .
A generic potassium current, I K , meant to represent an amalgam of the full family of potassium currents, is taken to have the form with E K = −90 mV (other parameters listed in Table A1).
The model could be further modified to include various Na + currents, such as a TTX independent sodium leak (Khaliq and Bean, 2010). Here, we included a TTX dependent persistent sodium current (Khaliq and Bean, 2010). Though not strictly necessary, its presence simplifies parameter selection for the activation variables. The persistent sodium current is taken to be of the form: where g pers = 0.002. Additionally, we include a general leak current of the form: I leak = g L (E leak − V) with E leak = −50 mV.

Calcium and Calcium-dependent Dynamics
The L-type calcium current is taken to have the form with E Ca = 100 mV and G L Ca given by, The calcium dependent SK potassium current is given by where G SK (u) =ḡ SK u 4 u 4 + K 4 1 . Calcium enters the cell predominantly via the L-type calcium channel (with a minor contribution due to the NMDA channel). The Ca 2+ is then taken to be ejected via a pump (akin to Canavier et al., 2007). The calcium, u, then varies according to where the cytosolic buffering constant is given by f Ca = 0.01; H is a lumped term involving the valence of calcium and Faraday's number (H ≈ 0.0193); the radius, r, of the neuron is taken to be r = 20 µm, and the pump dynamics are taken to be: with M pump = 500 (nM µm)/ms and K pump = 500 nM. Note that since intracellular Ca 2+ buffering is strong, u can be considered a slow variable.
1.1.4. Synaptic Currents: GABA, AMPA, and NMDA The two principal regions of the midbrain associated with DAergic neurons, the VTA and the SNc, both receive substantial amounts of GABAergic inhibition (from the ventral tegmental tail and the substantia nigra pars reticulata, respectively). The GABA cell populations display high frequency, tonic firing patterns and serve to hyperpolarize the target DA cell membrane potential. Since these rates are high, we consider an average conductance parameter for the inhibitory drive, that is, the fine-scale, detailed firing pattern of the GABAergic input is not explicitly tracked, and the current is taken to be: where g GABA is the conductance and E GABA = −65 mV is the reversal potential. In effect, this current acts akin to a very strong leak current. The DA cell receives excitatory, glutamatergic drive mediated by AMPA and NMDA receptors. The received AMPA current is taken to be I AMPA = g AMPA (E AMPA − V), where the reversal potential is E AMPA = 0 and the conductance evolves according tō where c = 0.002, σ s is spike strength, and χ noise is a binary parameter that determines whether noise is included. The spike times t i are generated via a Poisson process with a rate λ (the frequency for the noise typically taken between 10 and 50 Hz).
Computational studies with the inclusion of noise produce firing patterns that are more comparable to in vivo data that is inherently noisy.
The synaptic influence due to an incoming spike at time t i is subject to the α function: with the time constant given by τ α = 4 ms and represents a heaviside function. The NMDA current, also, influences DA neuronal dynamics. The NMDA current takes the form: with where Mg denotes the amount of magnesium, taken to be 0.5 µM.
Here E NMDA = 0,ḡ stim NMDA = 0.1,ḡ c NMDA = 0.01, and m e = 0.08. For simplicity, the parameterḡ stim NMDA is set to zero, except during NMDA activation/application. The form of the NMDA conductance depends upon the level of magnesium, shown below in Figure 3, where an increase in Mg 2+ shifts the conductance curve to be active at higher voltages.

Model Results
In the previous section, we presented a pared down representation of a DA neuron that includes an array of observed currents. In the following section, we demonstrate the model's ability to capture characteristic dynamics of DA neurons. By doing so, we identify conditions for eliciting a variety of firing modes and make a number of experimental predictions. In addition, we examine conditions akin to NMDA bath application (Johnson et al., 1992) for inducing burst firing and conditions akin to DA cell disinhibition (i.e., a release of inhibition from the GABAergic neurons).

Pacemaking and Burst Firing
Experimentally (Ping and Shepard, 1996) observed that application of the drug apamin induces burst firing in slice preparations by blocking the SK-channels, which reduces the flow of K + , that is, DA cell dynamics depend critically on the strength of the SK-conductance. Motivated by these results, we simulated this experiment by varying the strength of the parameter χ APA . For χ APA = 1, this corresponds to a strong SK-conductance and the DA cell fires tonically ( Figure 4A). However, as χ APA decreases, the neuron enters a weak SKconductance regime and the firing pattern transitions from tonic firing to endogenous or "natural" bursting ( Figure 4B). An initial tight spike doublet occurs at burst onset that is characteristic of DA cell activity (Grace and Bunney, 1984a). As expected, the strength of the SK conductance strongly influences a DA neuron's burst firing characteristics. , the neuron exhibits a tonic firing pattern (χ APA = 1). In slice preparations, application of the drug apamin suppresses the SK-conductance and yields an increase in burst firing. The mathematical model responds in a comparable manner, by reducing the SKconductance (χ APA = 0.2), the DA neural activity pattern enters a burst firing regime (a representative voltage trace is shown in B). In (C,D), we consider the calcium dependency behind these firing patters by constructing bifurcation diagrams for strong or weak SK conductances, shown in (C,D), respectively. We treat the intracellular calcium level as the bifurcation parameter, and as intracellular calcium, u, varies from high to low, the hyperpolarized stable steady-state (thick solid line) is lost via a fold bifurcation (with the unstable steady-state branch shown with a dashed line). For low calcium levels, the dynamics enter a periodic orbit and the neurons exhibit fast firing. The extrema values of the oscillations are represented by circles. Trajectories of the full dynamics are superimposed. Note that I 0 = 0 and that the dynamics for a weak SK conductance (χ APA = 0.2) admit an additional spike during a burst compared to when I 0 = 0.2, as shown in (B).
As the strength of the SK conductance affects DA cell firing patterns and the SK conductance is dependent on intracellular calcium levels, DA neural firing patterns are in turn dependent on intracellular calcium levels. The intracellular Ca 2+ buffering is strong, rendering u a slow variable. We treat intracellular calcium, u, as a bifurcation parameter in order understand the firing dynamics over a range of Ca 2+ levels, see Figures 4C,D.
If we consider ǫ = 2f Ca r ≪ 1 to be a small parameter where f Ca represents the cytosolic buffering constant and r represents the radius of the neuron. As ǫ → 0, the intracellular calcium level becomes approximately fixed. For "large" levels of calcium, the membrane potential attains a stable hyperpolarized steady state. As u decreases, the stable steady-state is lost through a fold bifurcation. For "low" Ca 2+ levels, the neuron enters a fast firing regime. This makes intuitive sense as for low Ca 2+ , the SK-channel remains inactive, while for high Ca 2+ , the SKchannels remain open resulting in a hyperpolarized state due to potassium leaving the cell in a non-gated fashion. When the full dynamics are included, starting from a hyperpolarized state, the neuron slowly becomes depolarized; simultaneously, the intracellular calcium level slowly decreases. Eventually the trajectory, passes the fold bifurcation and enters a periodic orbit associated with neural firing. In this depolarized periodic orbit, intracellular calcium levels rises. When intracellular calcium levels rises sufficiently to activate the SK potassium current, the trajectory falls off of the periodic orbit and returns to a hyperpolarized "resting" state. The activation of the SK-current coincides with the periodic branch colliding with the unstable steady state branch at a saddle-node on a limit cycle (SNIC) bifurcation. The shape of the bifurcation diagram is dependent on the relative strength of the SK-conductance (χ APA ). For strong SK-conductances, the neuron fires once (twice in some instances) before returning to a rest state, whereas for weak SK-conductance, the "bend" in the bifurcation diagram is more pronounced and there is a wider range of calcium levels where the periodic orbit occurs before colliding with the steady state branch. As a result, for weak SK conductance, the neuron goes through a greater number of oscillations in that depolarized oscillatory orbit before falling off the orbit and returning to the hyperpolarized state, i.e, the interburst quiescence. In addition, the interspike intervals (ISIs) increase through burst duration as is indicative of classical square-wave-like dynamics (Izhikevich, 2000). In future studies, we will be further examining the structure of the dynamics by performing fast/slow analyses with averaging.

Depolarization Block
The introduced model for DA cell activity undergoes depolarization block when overly driven, and we construct bifurcation diagrams to study the role of the amalgam term, I 0 , that represents contributions due to an applied current or bulk drive to the neuron, see Figure 5. For small values of I 0 , the neuron (for both weak and strong SK conductances) has a stable hyperpolarized steady-state, and the neuron exhibits quiescence. As I 0 increases, the steady-state loses stability through either a subcritical Hopf or through a fold bifurcation for strong and weak SK-conductances, χ APA = 1, 0.2, respectively. As I 0 increases, the frequency of periodic orbits increases, that is, the neuron fires with an increasing frequency. For sufficiently large I 0 , the periodic orbit is lost and the neuron enters as stable steady-state, associated with depolarization block.
For weak SK conductance (χ APA = 0.2) and the same progression from low to high levels of I 0 , the neuron begins at a hyperpolarized steady-state and enters a burst firing mode, that is a mixed mode oscillation at I 0 ≈ 0. As I 0 further increases, the frequency of the burst firing increases. At I 0 ≈ 1, the burst firing pattern transitions to a fast firing tonic pattern, where the mixed-mode oscillations appear to be lost. Past I 0 ≈ 3.5, neurons exhibit depolarization block and attain a stable depolarized steady-state solution. The exact nature of the substructure of the periodic firing patterns displayed for weak SK neurons and the transitions between the mixed-mode to fast firing is not examined in this work but is to be the subject of a future study. Beyond exhibiting multiple activity modes, a comparison of the bifurcation diagrams for strong and weak SK conductances shows that neurons with SK conductance have periodic orbits (i.e., fire) for a wider range of I 0 . This suggests that neurons with weak SK-conductances are more resistant to depolarization block, that is, neurons with a weak SK conductance will require higher levels of I 0 before a cessation of firing. Moreover, neurons with an extremely weak SK conductance and receiving high levels of drive may enter a tonic fast firing regime. A strong SK conductance appears to be the key to property of depolarization block observed for many DA neurons.

Burst Termination via Depolarization Block
The model exhibits a wide range of firing modes/behavior dependent on the make-up of an individual neuron's ionic currents, to say nothing of the effects of differences in broader system level activity changes in the neural circuit. For DA cells, classic burst firing is initiated by a spike doublet followed by intraburst spikes with ISIs that are approximately equal as classified in Grace and Bunney (1984b). However, burst firing patterns of a markedly different structure have been observed, e.g., even in the same (Grace and Bunney, 1984b) paper. Some of the burst firing patterns display decreasing spike amplitudes and increase of firing frequency during the burst event, so that the ISI decreased within a burst. These burst events appear to terminated via depolarization block. This type of bursting pattern was also found in the classic paper of Ping and Shepard (1996) in SNc DA neurons in slice preparations under apamin application and has been viewed in VTA extracellular recordings (Brischoux et al., 2009;Exley et al., 2011).
It is regarded as a characteristic of DA neurons that they are susceptible to depolarization block, whereby an increase of an applied current results in a cessation of DA cell firing due to sodium inactivation. In this work, we have identified an alternate set of parameters where individual neurons exhibit depolarization block in their endogenous dynamics, that is, the bursts terminate via depolarization block. The changes to the parameter regime include a broadening of the region of sodium activation, strengthening and altering the delayedrectifier potassium current, and changing the dynamics of the sodium inactivation gating variable. The detailed changes of parameters are listed in the Appendix (Table A2).
With these parameter changes to our model, the observed dynamics display the following described behaviors. For moderate levels of I 0 , neurons with a strong SK conductance exhibit regular firing at relatively high frequencies for DA neurons. However, for the same level of I 0 , weak SK-conductance neurons exhibit endogenous burst firing where the bursts are A B FIGURE 5 | Bifurcation diagrams for understanding how DA firing patterns depend upon the variable I 0 for strong or weak SK conductances, (that is, χ APA = 1.0 and χ APA = 0.2), shown in (A,B), respectively. Here we plot membrane potential with the generalized intrinsic/applied current variable I 0 taken as the bifurcation parameter. In both cases (weak and strong SK conductance), for small I 0 , the neuron exhibits quiescence. For increased I 0 , the dynamics follow a periodic orbit, i.e., the neuron fires, repetitively. For larger I 0 , the neuron undergoes depolarization block. The level of I 0 required to reach depolarization block depends inversely on the strength of the SK conductance. Note that for weak SK conductance for I 0 between approximately 0 and 1.0 that the periodic orbit displays a mixed-mode oscillation, indicative of endogenous burst firing, as I 0 increases, the firing displays "fast" pacemaker firing. Stable steady states are shown by thick lines and periodic solutions are labeled. terminated via depolarization block, shown in Figures 6A,B. Moreover, both the spike amplitude and ISI decreases throughout the burst duration, and burst termination coincides with a prolonged depolarization plateau. We plot these trajectories in phase space with respect to the membrane potential, V, and sodium inactivation, h, see Figures 6C,D. We note that for a weak SK conductance (χ APA = 0.2), a neuron, initially at a hyper polarized state (labeled H), becomes slowly depolarized and then enters a burst firing event with 4 spikes. During this burst event, both the spike amplitudes and the ISI decrease, that is, the spikes get smaller and there is an increase in the frequency during the burst. At the end of the burst event, the trajectory is at the point labeled D. Here the sodium inactivation is small and the trajectory lingers at the intersection of the V and h nullclines, an unstable steady state associated with long the plateaus of depolarization.
In this modified parameter regime, Ca 2+ buffering is fast, rendering Ca 2+ dynamics slow. We take the intracellular calcium to be at quasi-steady state and generate a bifurcation diagram in order to better understand the dynamics of the full system, shown in Figures 6E,F. The bifurcation structure contains a fold and a subcritical Hopf. For the case of a strong SK conductance, a hysteresis loop for regular firing exits whereby a single spike occurs followed by entering the hyperpolarized state. As Ca 2+ decreases, the trajectory falls from the fold and another spike occurs. For low levels of Ca 2+ , there exists a stable depolarized steady-state. In the case of a weak SK conductance, the shape of the bend in the curve is more pronounced and there is a larger regime for the periodic branch (this is an oscillation in voltage and corresponds to neural firing). As the neuron fires, Ca 2+ ramps up in the cell leading to the activation of the SK current. As this occurs on the bifurcation diagram, the trajectory would be approaching the subcritical Hopf bifurcation from left to right. When passing this bifurcation, the trajectory lingers on the unstable (depolarized) fixed point before dropping to the lower stable branch (i.e., the hyperpolarized state).
As the neuron loses Ca 2+ , the neuron becomes slowly depolarized. Eventually, it loses the stable steady-state via a fold bifurcation and enters the periodic orbit again, associated with neural spiking. During this spiking stage, intracellular Ca 2+ rises and the trajectory moves to the right in the bifurcation diagram. As this occurs, the amplitude of the spikes decreases during the burst progression, as also do the ISIs (that is, the spike frequency increases during a burst event). At burst termination, the trajectory lingers at the unstable depolarized steady state before "falling" down to the hyperpolarized interburst steady-state. Another difference for this parameter choice is that both the strong and weak SK conductance neurons have similar depolarization block properties with respect to the parameter I 0 , as seen in Figure 7. Both strong and weak SK conductance neurons are more impervious to increases in I 0 .

Quantifying Bursts
The seminal work of Grace and Bunney (1984a) examined the properties of in vivo bursting and introduced criteria for identifying the bursts. They noted that burst onset was associated with a tight doublet of two successive spikes, and thus identified burst onset by "the concurrence of two spikes with an interspike interval (ISI) of less than 80 ms. Burst termination was defined as an ISI of over 160 ms." This is a straightforward, algorithmic definition of a burst firing event and is still widely utilized. With this definition, one often employs a measure of the degree of burst firing that corresponds to the percentage of spikes occurring within bursts (spikes within bursts, SWB). Due to its prevalence in the literature, we consider this familiar burst measure (SWB), however with this definition of bursting with increased neuronal activity, it classifies fast, tonic firing to be equivalent to a burst of an extremely long duration. Thus, we contrast the SWB measure with an alternate burst measure introduced in van Elburg and van Ooyen (2004).
In order to differentiate, high rate, tonic firing from truly "bursty" dynamics, we consider the Van Elburg/Van Ooyen measure for degree of bursting. This measure utilizes both the interspike and "two-spike" interval. Take ISI to be the distribution of the ISIs over an experiment. In addition, let TSI represent the "two-spike interval" distribution for the same spike train, which is made up of time intervals corresponding to the time between a spike and the second subsequent spike. If σ I and σ T represent the standard deviation of the ISI and TSI distributions, respectively, and µ I is the average value of the ISI distribution. The Van Elburg/Van Ooyen measure for the degree of bursting is given by Using this measure, B above a threshold value, B , would constitute a bursting neuron. For example, in van Elburg and van Ooyen (2004), a neuron with a measure B > B = 0.15 was considered to be bursting.

DA Activity Modes
A wide range of parameter values were utilized in order to examine the necessary conditions to elicit certain modes of activity that have been observed by Hyland et al. (2002) and Mameli-Engvall et al. (2006) and others. Other laboratories are working to identify (and classify) DA neurons from electrophysiological properties (Shepard and German, 1988;Lammel et al., 2006;Margolis et al., 2006;Li et al., 2012). We examine the mechanisms behind this ensemble of dynamical classes of DA neurons. In Figures 1A,B, we presented the experimental ensembles of neural activity found by Mameli-Engvall et al. (2006) and compared with the firing modes generated from our computational model amount of drive and containing neurons with weak to strong SK conductance. The model results expose a high correlation of weak SK conductance with the high firing regime (both bursting and tonic). Moreover, neurons with weak to medium SK conductances are also associated with a significant portion of the burst firing events. In fact, it is the balance of L-type Ca 2+ to the SK Ca 2+ -dependent potassium current that shapes the activity profile. Throughout this paper, we consider a fixed level of L-type Ca 2+ and vary the SK conductance, however it would be equally valid to keep the SK conductance fixed and vary the L-type Ca 2+ . In addition to these two currents, the level of GABAergic inhibition plays a strong role, and we go on to study how it can suppress or sculpt the firing dynamics of DA neurons and their displayed activity modes.
Here we consider a range of strengths of the SK conductance and a range of values for the level of the intrinsic/background drive variable I 0 , shown in Figure 8, in order to study the mechanisms underlying multiple activity modes. From the phase diagram, we deduce that for moderate to strong SK conductance strengths, the DA neurons undergo depolarization block when given a large applied current, I 0 . The model results suggest the experimental prediction that neurons that have undergone apamin exposure or that have an inherently lower SK conductance will be more resistant to depolarization block (that is, require a greater applied current to initiate it). There appears to be an inverse relationship between SK-conductance strength and the needed applied current to bring on depolarization block. In addition, neurons that display higher firing rates appear to be correlates with lower SK levels.
By considering the two measures of bursting (SWB measure and the Van Elburg/Van Ooyen measure), underlying differences in the firing patterns can become apparent. The SWB measure classifies neurons with a weak SK conductance that fire tonically at approximately 8 Hz to be bursting, whereas the burst measure of Van Elburg and Van Ooyen does not classify these neurons as bursting. Moreover, the Van Elburg and Van Ooyen measure accentuates the presence of mixed-mode oscillations within the firing patterns that are associated with burst firing. Another difference between the measures is that at a level of I 0 where neurons are undergoing depolarization block, the neurons fire sporadically due to noise. At these points, the burst measure classifies these neurons to be bursting.
The strength of the SK conductance and the intrinsic/applied current do not alone determine a DA neuron's activity mode. The level of GABAergic drive plays an import role in determining the activity profile of a DA neuron as noted in Komendantov et al. (2004) and Lobb et al. (2010Lobb et al. ( , 2011. We generate a series of diagrams akin to Mameli-Engvall et al. (2006) to study the impact of GABA on influencing DA activity modes (we consider low, medium and high levels of GABA), shown in Figure 9. The most telling conclusion that we infer from these plots is that the strength of the GABAergic tone plays a dominant role in DA cell activity mode determination. For a strong GABAergic tone, burst firing is reduced and the activity modes are more tonic in nature. As GABAergic tone decreases, firing rates increase as well as the number of neurons exhibiting burst firing. This suggests that if the nicotinic acetylcholine receptors (nAChRs) within the GABAergic population are exposed to nicotine and enter a desensitized state, the resulting disinhibition of the DA cell population will result in a bulk increase of the average DA cell activity, but also an increase in burst firing. This increase in bursting occurs predominately in neurons that have weak SK conductances, so this supports the hypothesis of Mameli-Engvall et al. (2006) that neurons that display some degree of bursting are neurons that display increased burst responses to nicotine simulation, that is, there are neurons that are primed to burst.  (2004). For small I 0 the system displays quiescence. For "moderate" I 0 , DA cells display activity (neurons with a lower SK conductance fire at a moderately higher rate). For larger I 0 , the DA cells cease firing (that is they undergo depolarization block) and the level of I 0 required to initiate depolarization block is inversely proportional to the SK-conductance. That is, neurons with a low SK-conductance continue to fire at levels of I 0 that would cause a neuron with a high SK conductance to cease firing.
One matter to note is that with either measure of bursting, spike doublets are classified as bursts and much of the "low firing, high burst" activity are, in fact, doublets generated periodically. Differences in the degree of GABAergic control to a DA neuron greatly influences the resulting firing pattern. In rats under anesthesia via chloral hydrate, only 18% of the observed DA neurons in the SNc exhibited bursting, whereas, 73% of the VTA DA neurons exhibit burst firing (Grenhoff et al., 1988). In our numerical studies, we find that neurons receiving more GABAergic inhibition exhibit more pacemaker-like firing, which is more indicative of SNc DA neurons. Thus, we postulate that the relative degree of GABA inhibition to the SNc DAergic neurons is greater than that to the VTA DAergic neurons.

A Possible Novel Firing Mode Unveiled
In our numerical studies, numerous phase diagrams were generated considering many sets of parameters and amounts of noise in the glutamatergic synaptic drive. For a wide range of parameters, phase diagrams akin to the one pictured in Figure 9 were obtained. However, we also discovered another parameter set (listed in Table A3) that exhibited a number interesting characteristics. The associated phase plots for this parameter set are shown in Figure 10, where the activity modes for a number of FIGURE 9 | The strength of the GABA inhibition has a strong influence on the firing and burst firing characteristics of a DA neuron over a range of neuronal drive levels and SK-conductance strengths. Neurons with weak, moderate, and strong SK conductance are denoted via the use of squares, disks, and asterisks, respectively. Here GABA levels listed as low, mid, and high corresponds to g GABA = 0.01, 0.02, and 0.03, respectively. The GABA conductance included is taken to be constant and does not capture the dynamic network level changes in the mean field activity of the GABA neuronal population innervating the DA neuron.
GABA levels are considered as well as a range of SK conductances and values of I 0 .
We observe that there is a pronounced sensitivity to the level of GABAergic inhibition, whereby for high GABA levels, much of the neuronal activity is silenced with the burst firing modes almost completely removed. Furthermore, the threshold for I 0 to induce depolarization block varies greater at the different levels of GABA inhibition. Akin to the previous phase diagram, the firing rate of the neurons increases with a decrease of the SK conductance. However, the "very high" firing rates plotted here not only rise above 10 Hz, but as I 0 raises past 3, the firing rates enter fast firing regimes between 20 and 40 Hz. Such a tonic, fast firing regime is interesting for a number of reasons. For moderate levels of GABAergic inhibition (Mid GABA), for neurons with low SK conductance (0.15 < χ APA < 0.2), as I 0 increases, these neurons cease firing but with a further increase of I 0 , these neurons enter a fast firing regime. This fast firing regime may occur for extreme parameter values, but it may signify a subclass of DAeric neurons that defy the convention that DAergic firing rates are less than 12 Hz, as has been suggested by the findings of Lammel et al. (2006). These computational findings taken together with the data of Lammel et al. (2006) suggest the possible existence of a novel activity mode yet to be unveiled and yields to the experimental prediction that neurons exposed to apamin may display similar physiological dynamics. In Figure 11, we go onto examine the mathematical structure of this novel firing mode. We observe the emergence of a periodic solution, i.e., very fast neural firing.
FIGURE 10 | An alternate phase diagram that suggests a novel activity mode of DAergic neurons. The strength of the GABA inhibition has a strong influence on the burst firing characteristics of a DA neuron over a range of neuronal drive levels and SK-conductance strengths.
Neurons with weak, moderate, and strong SK conductance are denoted via the use of squares, disks, and asterisks, respectively. Again, GABA levels listed as low, mid, and high corresponds to g GABA = 0.01, 0.02, and 0.03, respectively.
FIGURE 11 | Bifurcation diagram for moderate GABA level with a weak SK conductance (χ APA = 0.2) to study the neural dynamics over a range of I 0 . The observed results are not in perfect accord to the phase diagram shown in Figure 10 for the frequency distribution at moderate GABA level as there is no inclusion of noise in this diagram.
However, we observe that as I 0 increases from 0, the neuron enters a firing regime, but as I 0 increases further it enters a state of quiescence. With a further increase of I 0 the neuron can enter another fast firing regime. Stable steady states are shown by thick lines and periodic solutions are labeled.

NMDA-induced Firing
Other than endogenous burst firing behavior, increases in the NMDA conductance elicits transient or endogenous burst firing. This is inline with the observed endogenous bursting patterns in vitro (Johnson et al., 1992) due to NMDA-bath application and is consistent with the NMDA-initiated burst firing in vivo (Chergui et al., 1993). Li et al. (1996) first modeled NMDA-induced bursting using a two-compartmental representation of the the neuron, which relied on somatodendritic interactions in order to generate burst firing. Here, we utilize the original parameter set in our single-compartment model (with sodium activation slightly altered, p 3 = 12.5, though these properties are general over a range of different parameters). We consider a range of SK conductances and different levels of Mg 2+ (different levels of Mg 2+ alter the activation dynamics of the NMDA current, see Figure 3). In Figure 12A, a DA neuron with a "strong" SK conductance and relatively low Mg 2+ level, transitions from a tonic firing regime to a transient burst of activity at time 2.0 s by an instantaneous increase of the NMDA conductance, g NMDA , by 0.1 ms/cm 2 for a duration of 2.0 s. After the burst coincides with the onset of the stimulus (sudden rise of the NMDA conductance) and then enters a depolarized state. (Note that whether the neuron enters a fully depolarized state depends upon the variables associated with sodium activation and inactivation). At the end of the stimulus, the neuron becomes hyperpolarized and displays a long recovery period before returning to pacemaker-like firing. For the same level of Mg 2+ (0.5 µM), a DA neuron with a weaker SK conductance can be similarly be induced to burst via an increase of the NMDA conductance (see Figure 12C). The resulting burst firing event is of a longer duration (that is, "burst firing is augmented, not suppressed, by SK channel blockade" Deister et al., 2009) and here is followed by pacemaker-like firing for the duration of the increased NMDA level.
For a similar rise in the NMDA conductance but with high levels of Mg 2+ (3.2 µM as in Johnson et al., 1992), at stimulus onset an initial burst event is induced (whose duration depends on the strength of the SK conductance). Following that initial burst, the neuron enters a burst firing regime (the degree of bursting dependent on the strength of the SK conductance) as shown in Figures 12B,D. In effect, with high levels of Mg 2+ , the NMDA conductance varies more during the observed rhythm. When the neuron is in a "down" state, the NMDA conductance is small and this accentuates the rhythm in part due to the interaction between the L-type Ca 2+ and the SK-current. Moreover, during the plateaus of depolarization the NMDA A B

C D
FIGURE 12 | Differential burst firing events due to an increase in the NMDA conductance for I 0 = 0.3, responses depend on both the strength of the SK conductance and the Mg 2+ level. In (A), a neuron with a relatively strong SK conductance (χ APA = 1.0) has NMDA applied at time t = 2 s (manifested by an immediate rise of the NMDA conductance), which elicits a burst firing event. In (C), the preparation is repeated but with a weak SK conductance (χ APA = 0.2), results in a significantly longer burst firing event followed by pacemaker-like firing for the duration of the stimulus. In (A,C) the Mg 2+ level was 0.5 µM, in (B,D), we repeat the experiment at an elevated level of Mg 2+ (3.2 µM) and note the resultant endogenous burst firing during stimulation (akin to Johnson et al., 1992). In addition to tracking the membrane potential, we also display the intracellular calcium dynamics. These simulations include noise in the AMPA drive of 40 Hz.

A B C D
FIGURE 13 | Burst firing events due to NMDA increases for neurons of strong and weak SK conductance (χ APA = 1.0 and 0.2) shown in (A,C) for low Mg 2+ levels, respectively from depolarized (or near depolarized) states (driven by I 0 = 2.5). NMDA increases by 0.1 at time t = 2 s for a duration of 2 s. In (A), there is a transient 4 spike burst, while for a weak SK conductance (shown in C) we observe a prolonged 9 spike burst event followed by sporadic activity until the NMDA stimulus is released. In both instances, at the release of stimulation, there is a pronounced hyperpolarization of the membrane potential. In (B,D), we consider similar conditions but for high levels of Mg 2+ . We find that in each case (for both weak and strong SK conductance) that the rise in the NMDA conductance at time t = 2 s results in quasi-regular firing.
conductance is active. We note that in the low Mg 2+ regime (shown in Figures 12A,C) that the NMDA current is present and driving the neuron, and with the rise in the conductance level, this NMDA current acts akin to a rise in I 0 . Please note that even for neurons receiving high levels of GABA with completely suppressed activities that NMDA can still initiate transient burst events (figure not shown).
Recall that the phase diagrams suggest that SK-weak DA neurons require higher levels of I 0 before they undergo depolarization block. These results suggests that the neurophysiological response to an NMDA application can be used to differentiate weak vs. strong SK DAergic neurons, that is, the relative strength of the SK conductance could be elucidated physiologically via differential excitable responses. Our finding suggests that in addition to weak SK neurons displaying a wide range of endogenous burst firing regimes, they also respond to stimuli in a more "bursty" fashion. This is line with the observation of Mameli-Engvall et al. (2006), that in response to nicotine, initially bursting neurons are more apt to enter stronger burst firing regimes, that is, it is unlikely that initially tonic firing neurons will enter burst firing regimes.
We consider the neural response to bulk rises in the NMDA conductance by DA neurons receiving a strong applied current (that is, we consider a relatively large value for the parameter I 0 = 2.5). For both low and high Mg 2+ , when the drive is high I 0 , neurons with a strong SK conductance are near depolarization block (with sporadic firing due to noise) and neurons with a weak SK conductance exhibit "fast" pacemaker firing, see Figure 13. For low Mg 2+ (Figures 13A,C), an instantaneous rise in the NMDA conductance at t = 2 s elicits a similar response to the experiment performed with moderate I 0 (I 0 = 0.3). Akin to the previous simulation at moderate levels of I 0 , at the completion of the stimulus the cell enters a hyperpolarized state. By and large, there is a decrease of activity due to the prolonged raise in the NMDA levels. For higher levels of Mg 2+ , during the period of elevated NMDA, the neurons display fast firing (that is, elevated NMDA elicits an initial burst then leads to an increased firing rate for both weak and strong SK conductances, Figures 13B,D. For high levels of Mg 2+ , the degree of hyperpolarization following the NMDA stimulus is diminished.

Disinhibition-Induced Firing
Through an examination of the phase plots of the activity modes for various levels of GABA, we noted the profound effect of GABA on influencing the firing modes. However, that study was for the endogenous firing patterns of neural dynamics. Here we consider transient disinhibition events (akin to Lobb et al., 2011). From a state that exhibits low neural firing due to the degree of high GABAergic tone, a release from inhibition results in transient burst firing event followed by repetitive firing in a slow tonic or burst firing pattern, dependent on the the SK-conductance of the neuron, Figure 14. From an initial state of high GABA inhibition, g GABA = 0.04, neurons are released from inhibition at t = 3 s (g GABA = 0.04 → 0.01), inducing transient burst firing events at the onset of disinhibition (whose duration depends on the SK conductance and the level of Mg 2+ ). In (A,C), the Mg 2+ level is 0.5 µM, and there is some background activity both prior to and post disinhibition events. In (A), a neuron with a relatively strong SK conductance (χ APA = 1.0) displays a brief transient burst at the onset of disinhibition then enters a pacemaker-like firing regime. Whereas in (C), a neuron with a weak SK conductance (χ APA = 0.2) yields a longer burst firing event at the onset of disinhibition followed by endogenous burst-like firing. In (B,D), we repeat the experiment at an elevated Mg 2+ level and consider high and low SK conductances (χ APA = 1.0, 0.2), respectively. In (B,D), there is nearly no background activity prior to and post disinhibition events, and the resulting burst firing events at disinhibition onset are significantly longer than the case for low Mg 2+ levels, shown in (A,C). In (A-D), following the period of disinhibition, the cell returns to its previous firing pattern with no appreciable post-stimulus hyperpolarization.
The simulated disinhibition event was represented via changes in the conductance level for the GABA, where high GABA levels transition to low GABA during the disinhibition event, that is, g GABA = 0.04 → 0.01) during the disinhibition event, that is, this disinhibition represents a decrease in the shunting current, rather than an increase in a constant drive current, say I 0 . The model results, albeit general in line with, do show some differences in from Lobb et al. (2011). During disinhibition (relatively short episodes) in the former study, the neurons displayed increased fast firing, i.e., a burst events. In our studies, we consider longer periods of disinhibition than in the former study. At disinhibition onset, neurons respond with a burst firing event followed by a rhythmic firing pattern (which may display of degree of bursting dependent on the SK conductance of the neuron, see Figure 14C). For the parameters considered here, for high Mg 2+ levels with strong inhibition was near quiescent, so that neural activity was relegated to the period of disinhibition (see Figures 14B,D) as the NMDA current would contribute very little with Mg 2+ = 3.2 µM and V ≈ −50 mV (please refer to Figure 3). In addition, we note that after a disinhibition event there is no resulting hyperpolarization as found after an NMDA bath application. We also consider how the degree of disinhibition can influence the neural response at disinhibition onset, that is, the number of spikes. Our findings, displayed in Table 1 suggest a graded response with respect to the degree of the disinhibition event and the number of spikes elicited. The degree of disinhibition may represent detailed changes in the reward expectation, our results suggest a continuous, rather than all-or-none, response to disinhibition events. However, we note that for high SK conductance (i.e., χ APA = 1) and with a very weak degree of disinhibition, the elicited response is minimal.

Discussion
Here in this paper, we have utilized computational modeling in order to better understand the mechanisms controlling the ensemble of DA neuronal classes found in DA neurons: both endogenous rhythms and stimulus-induced firing patterns. These findings are in line with experimental findings and complement previous theoretical modeling efforts, however in a single compartment framework as somatodendritic interactions appear to not be a necessary for burst firing initiation (Blythe et al., 2009;Deister et al., 2009) (in particular the distal dendrites). By performing numerical simulations for a variety of parameter sets, our model leads to a number of experimental predictions. Our studies suggest that neurons with a relatively low SK conductance fire in a high firing mode, are correlated with burst firing, and require higher levels of applied current before undergoing depolarization block. We considered the role of GABAergic inhibition on an ensemble of neural types and examine alternate scenarios for both NMDA-initiated and disinhibition-initiated burst firing. Moreover, we may have unveiled another class of DA neurons that enters very fast firing regimes (though such a parameter regime may prove to be unphysiological).
This work accentuates the importance and roles of the SK conductance, GABA inhibition, and the NMDA current have in burst generation. Our results suggest that neurons with a relatively low SK conductance fire in a fast firing mode, are correlated with burst firing, and require higher levels of applied current before undergoing depolarization block. Due to the interplay between the L-type Ca 2+ and the SK conductance, another way to consider these low SK conductance neurons would be strong L-type Ca 2+ neurons. However, a decrease of the SK conductance is of potential physiological significance since such slow K-channels are the primary target for downmodulation by ACh through muscarinic AChRs (mAChRs), hence ACh action through mAChRs may prime the DA neurons for bursting (e.g., Drion et al., 2010). Further importance of the potassium currents is evidenced by an alteration of the potassium dynamics, which we found can lead to neural bursting that terminates via depolarization block, akin to firing pattern observed in vivo in VTA DA neurons and endogenous burst dynamics in SNc DA cell preparations under apamin application. The spikes within these bursts are of decreasing amplitude and have decreasing ISIs. This suggests a markedly different bifurcation structure behind the dynamics. Due to the existence of two characteristic burst firing patterns, we postulate that this may be indicative of VTA DA bursts arising due to two different methods: NMDA driven vs. a release from GABA inhibition.
Our computational studies on NMDA and GABA are qualitatively in line with the findings of Komendantov et al. (2004) in the context of multi-compartmental models and with the study of Lobb et al. (2011). The firing patterns we observed were characterized by an initial burst event and then enters either a firing mode dependent on the characteristics of the cell (that is, not a sustained high frequency firing for the duration of the disinhibition). The model reinforces the framework that NMDA levels, strength of the SK-conductance, and the level of GABAergic inhibition are all important components in determining the firing mode of a DA neuron, however it further refines the roles of the individual components. The level of GABAergic inhibition strongly influence the endogenous burst firing properties of neurons, that is, strong GABA inhibition leads to mostly tonic firing activity modes. Whilst a removal of that inhibition, i.e., disinhibition, can induce transient burst disinhibition events. Moreover, the main initiator of burst firing events is likely NMDA. We can differentiate NMDA stimulation vs. disinhibition events by the subsequent drop in the membrane potential after prolonged NMDA stimulation. The activation properties of the NMDA receptors can dramatically affect the response properties of neurons. In our studies, we considered a range of NMDA conductances and vary the amount of Mg 2+ . Note that references to high or low amounts of Mg 2+ are equivalent to shifting the activation profile for the NMDA receptors (see Figure 3). In NMDA stimulation experiments, DA neural response was characterized by a transient burst event followed by either quiescence or pacemaker firing in the case of low Mg 2+ , whereas for high Mg 2+ , DA neurons displayed pronounced bust firing or high firing rates. The response properties of NMDA stimulation differed from disinhibition events whereby neurons displayed an initial burst at disinhibition onset followed by either pacemaker-like or a burst firing for neurons with strong and weak SK conductances, respectively. Taken together, the above results suggest to obtain strong endogenous bursting under NMDA stimulation, neurons should have both weak SK-conductances and high levels of Mg 2+ , whereas for disinhibition events, the strength of the SK conductance alone appears to account for endogenous burst firing patterns.
In considering the role of the SK conductance and GABAergic tone, we identified a set of parameters whereby neurons with weak SK conductance could undergo a cessation of neural firing via an increase in I 0 , however with further increases of I 0 would enter a fast spiking firing mode. This puts into question the ubiquitous belief that DA neurons fire at a rate less than 12 Hz. This suggests that here may exist a subclass of DA cells that exhibits fast firing as proposed in Lammel et al. (2006) where they the identified "a type of dopaminergic neuron within the mesocorticolimbic dopamine system with unconventional fastfiring properties." In Lammel et al. (2006), the observed fastfiring DA neurons were able to sustain higher firing rates above the conventional upper limit of 10 Hz. It would be interesting to have measurements of the differences of the relative strengths of the L-type Ca 2+ and SK conductances for the fast firing neurons observed by Lammel et al. Moreover, it would also be of interest to examine the neural firing patterns of VTA DA neurons from slice preparations with apamin administered. With increased applied current, would these neurons similarly or reliably enter a fast firing mode?
The framework we have put forth represents a step forward in capturing DA cell dynamics and could be further improved to include more detailed second messenger signaling processes in order to elucidate further differences between the two principal classes of midbrain DA neurons: those of the VTA and the SNc. However, the GABA studies in this paper, suggest that the SNc DA neurons likely receive a higher level of GABA inhibition due to the observation that SNc neurons display more pacemaker-like firing patterns than VTA DA neurons. The model framework itself could be extended from this conductance-based single compartment model of DA cell activity to include the detailed dynamics of two prevalent nAChR subtypes (α4β2 and α7) and the averaged-activity of the GABA cells (as developed by Graupner and Gutkin, 2009). This new framework would take into account a multitude of scales (i.e., linking the VTA mean field activity, detailed conductance-based representations of DA neurons, and the dynamics of the nAChRs) and could yield further insight to the mechanisms underlying burst firing in vivo and aid in elucidating the biophysical manifestations of nicotine addiction.

Author Contributions
The paper was developed by AO and BG with PF providing valuable biological insight. The paper was principally drafted by AO with input for revisions from both BG and PF. All numerical studies were performed by AO.