Interaction of NMDA Receptor and Pacemaking Mechanisms in the Midbrain Dopaminergic Neuron

Dopamine neurotransmission has been found to play a role in addictive behavior and is altered in psychiatric disorders. Dopaminergic (DA) neurons display two functionally distinct modes of electrophysiological activity: low- and high-frequency firing. A puzzling feature of the DA neuron is the following combination of its responses: N-methyl-D-aspartate receptor (NMDAR) activation evokes high-frequency firing, whereas other tonic excitatory stimuli (-amino-3-hydroxyl-5-methyl-4-isoxazolepropionate receptor (AMPAR) activation or applied depolarization) block firing instead. We suggest a new computational model that reproduces this combination of responses and explains recent experimental data. Namely, somatic NMDAR stimulation evokes high-frequency firing and is more effective than distal dendritic stimulation. We further reduce the model to a single compartment and analyze the mechanism of the distinct high-frequency response to NMDAR activation vs. other stimuli. Standard nullcline analysis shows that the mechanism is based on a decrease in the amplitude of calcium oscillations. The analysis confirms that the nonlinear voltage dependence provided by the magnesium block of the NMDAR determine its capacity to elevate the firing frequency. We further predict that the moderate slope of the voltage dependence plays the central role in the frequency elevation. Additionally, we suggest a repolarizing current that sustains calcium-independent firing or firing in the absence of calcium-dependent repolarizing currents. We predict that the ether–a-go-go current (ERG), which has been observed in the DA neuron, is the best fit for this critical role. We show that a calcium-dependent and a calcium-independent oscillatory mechanisms form a structure of interlocked negative feedback loops in the DA neuron. The structure connects research of DA neuron firing with circadian biology and determines common minimal models for investigation of robustness of oscillations, which is critical for normal function of both systems.


Introduction
The central dopamine system is involved with many behavioral and cognitive tasks [1]. Dopamine neurotransmission is found to play a role in addictive behavior [2] and is altered in psychiatric disorders [3]. The DA neuron fires a burst of activity whenever an animal receives an unpredicted reward or a reward-predicting stimulus [4]. The mechanisms underlying the excitability of the DA neuron have attracted the interest of electrophysiologists for many years. One of the major questions is what synaptic input triggers the excitation, or is a combination of inputs required? The stereotypic response of the DA neuron to the unpredicted reward is a burst of activity with firing frequencies above 20 Hz. In vivo, the bursts are singular events superimposed on a low-frequency single-spiking background firing (1)(2)(3)(4)(5). A burst is shown to be triggered by the stimulation of glutamatergic inputs, of which NMDA, but not (AMPA) receptor activation was found to be critical [5][6][7]. In vitro, these distinct responses to the activation of the two types of glutametergic synapses are confirmed: bursts evoked by stimulations in slices are dependent on the activation of NMDA receptors, whereas tonic activation of AMPA receptors alone does not elicit bursting [8,9]. This is a puzzling feature of the DA cells that distinguishes it from many other types of neurons: How does DA neurons differentiate the NMDAR and AMPAR stimulation?
Similar to the AMPAR stimulation, tonic somatic current injections cannot elevate the firing frequency significantly: Applied depolarization blocks firing in the DA neuron as the frequency increases above 10 Hz [10,11]. Interestingly, somewhat higher frequencies are possible during a transient after the onset of the current injection [10,11], but not after the onset of a tonic AMPAR stimulation [9], which evokes only a single spike. These low frequencies and the susceptibility to the blockade of firing are features that distinguish the DA neuron.
Our previous model [12] has suggested how the neuron can respond differentially to the applied current and AMPA vs. NMDA receptor stimulation. Several laboratories tested our major predictions, and found one of them inconsistent with the new experimental data [9,11,13]. Namely, the false prediction was the dendritic drive for the high-frequency firing: Because of the difference in morphology, the isolated dendrites must have a much higher natural frequency. We had found how parts of the neuron can dominate its dynamics and, consequently, determine the frequency. To change its frequency, the neuron needs only a stimulus that switches the dominance from one part of the neuron to another. We had shown that NMDAR switches the dominance from the low-frequency soma to high-frequency dendrites in the model. However, the switch in dominance required distribution of the NMDA stimulation onto fine dendrites, whereas somatic NMDA stimulation was unable to elevate the frequency. The new experiments have demonstrated the opposite -somatic NMDAR stimulation or a virtual current elicits high-frequency firing [9,13]. This disproves our mechanism.
In this paper, we design a new mechanism for the differential responses of the DA neuron to applied depolarization and AMPA vs. NMDA receptor activation. The major distinction with the old mechanism is that the gradient of the natural frequencies along the dendrites is not used. First, this allowed us to achieve the high frequencies either for somatic or dendritic, as well as whole neuron NMDAR stimulation. Consistent with experiments [11], a proximal or somatic stimulation is more efficient in the model. Second, it allowed for the major reduction of the model to only one compartment. Based on this reduction, we give a very simple explanation of the frequency elevation. Remarkably, the biophysical properties described by the model are exactly the same as in the previous one [12]. However, the calibration of the model is different, and the most important is an altered voltage dependence of the NMDAR current. Consistent with experiments [9,13], our study highlights that the voltage dependence provided by the magnesium block of the NMDAR current is critical for the elevation of the frequency. We further predict that, of the parameters of the NMDAR current, the moderate slope of its voltage dependence is central in the new mechanism. The results refocus the attention of the experimentalists from the search of new determinants of the burst firing to the well-known oscillatory mechanism in the light if its new interaction with the NMDAR current.
Another surprising feature of the DA neuron accentuated recently is its calcium-independent firing: the neurons continue robust firing at a low pacemaking rate after pharmacological blockade of Ca 2+ channels or Ca 2+ chelation, or in Ca 2+ -free solution [11,[14][15][16][17][18]. How can a neuron display calcium-independent firing if its oscillatory mechanism is consistently shown to rely on the calcium and calcium-dependent potassium currents (see e.g. [19][20][21])? This evoked a vivid discussion about the contribution of calcium vs. sodium channels to the depolarization phase of the voltage oscillations and resulted in the conclusion that both types of channels contribute to the same oscillatory mechanism [22] and that calcium influx is not essential for oscillations [13]. What receives no attention in these papers is that something must repolarize the membrane when the calcium-dependent potassium current does not activate. This cannot be a strong delayed rectifier potassium current because this current would support a high frequency indifferently to the type of stimulation: applied current, AMPA or NMDA receptor stimulation.
We suggest a voltage-dependent potassium channel that supports oscillations when the calcium-dependent repolarization mechanisms do not activate. No current had been suggested for this role, and we start with the assumption-free basic potassium current and calibrated it according to its function of supporting oscillations with particular characteristics. Then, we identify the current as the ether-a-go-go (ERG) current. ERG1 proteins are expressed in the substantia nigra pars compacta neurons [23]. Experiments show that the ERG current is functional in the DA neuron [24,25]. The new current allows us to reproduce another recent result -the high-frequency firing persists under an SK current blocker apamin [9]. We study the interaction between the repolarizing currents and suggest how calcium-dependent and calcium-independent oscillatory mechanisms are combined in the DA neuron. The currents form a structure in which two feedback loops share components, and we adopt the term of interlocked feedback loops from another discipline, circadian biology for our new representation of the DA neuron. This establishes a common framework for understanding robustness and regularity of oscillations critical in both disciplines.

Materials and Methods
We build a reduced model of the DA neuron to understand its properties rather than a detailed biophysical model. Thus, only currents that are shown to be critical for oscillations are included. The structure of the model is similar to that in our previous publication [12]. The major changes are altered NMDA voltage dependence and a newly introduced ERG-type voltage-dependent potassium current. The model is presented in two morphologies: (1) a reconstruction of a DA neuron and (2) a single compartment that ignores the spatial structure of the neuron. We show that these two morphologies display very similar patterns. The biophysical mechanisms included in both morphologies are the same. They are described below.

Biophysical Properties
Both one-compartment and reconstructed morphology models studied below implement the following membrane mechanisms. The central currents in the model are an L-type non-inactivating voltage-dependent calcium current, an SK-type Ca 2+ -dependent potassium current, and a voltage-dependent potassium current, calibration of which suggests that it's an ERG-type current. Gating of the calcium current is fast, and it is treated as instantaneous for simplicity. Gating of the Ca 2+ -dependent potassium current is also much faster than the timescale of its gating variable, Ca 2+ concentration. Thus, we omit an extra gating variable for this current, treating it as instantaneous as well. By contrast, gating of the ERG current requires an explicit gating variable for generating oscillations under blockade of the Ca 2+ -dependent current. We calibrate it to play this role and discuss its parameters as a model prediction.
The four resulting equations are: Here, v is voltage, ½Ca 2z is calcium concentration, n is the gating (activation) variable of the ERG current, and h is the inactivation gating variable for the fast sodium current. In the voltage equation (1), g Ca (v) is a voltage-dependent Ca 2+ conductance; g KCa (½Ca 2z ) is a Ca 2+ -dependent K + conductance; g ERG (n) is the ERG conductance. A small leak conductance g l and a voltage-gated instantaneous potassium conductance g K (v) are included to limit the input resistance and to confine the voltage in the subthreshold range, respectively. We include a spikeproducing fast sodium current with maximal conductance g Na . The current was found not to be essential for generating either high-or low-frequency activity [9,20]. Thus, we calibrated the model to display firing patterns similar to the subthreshold oscillations when the fast sodium current is blocked.
Gating for the fast sodium current is in the standard Hodgkin-Huxley form. The activation of the current is assumed instantaneous and described by the function.
The inactivation kinetics is described by Eq. (4), where.
The current is calibrated to produce a spike per each maximum of the subthreshold oscillations without qualitatively changing the oscillatory pattern.
The equation for the gating variable of the ERG current (3) is in the standard Hodgkin-Huxley form with the following activation function and time constant: The functions are shown in Fig. 1A, and their parameters are a product of calibration of this current to sustain pacemaking in the absence of the Ca 2+ -dependent potassium current. The properties of this pacemaking (see below) allow us to specifically determine the half-activation of this current, while the time constant is determined qualitatively rather than quantitatively. The parameter search was performed manually because the criteria are hard to automate.
Conductances of the currents in Eq. (1) are given by the following functions: , g KCa (u)~ g g KCa ½Ca 2z 4 ½Ca 2z 4 zk 4 , g KV (n)~ g g ERG n 4 , Here, all maximal conductances are given per unit of surface area in mS=cm 2 (see Table S1 for the values of the parameters). The voltage dependence of the Ca 2+ current replicates the characteristic low-threshold L-type current found in DA neurons [26,27]. The calcium current is shown to be non-inactivating [20]. The fourth power dependence on Ca 2+ concentration above is typically used to best represent the characteristics of the SK type Ca 2+ -dependent potassium current [28]. We choose its parameters in the range estimated for the DA neuron [20].
The additive parameter I in Eq. (1) represents an applied current (inmA). It is normalized by the surface area of the neuron A (in cm 2 ). The conductance g g AMPA is a constant conductance density of a linear AMPAR synaptic current. The nonlinear function of the voltage with g NMDA (v)~ g g NMDA 1z0:1½Mg 2z e {0:062v reflects the activation of the NMDAR synaptic current (see Fig. 1B). Here, ½Mg 2z is magnesium concentration, and the magnesium block is treated as instantaneous. The parameters of this voltage dependence are very close to the values previously used in the literature [12,29]. An important modification is a decrease in the slope of the NMDA conductance activation m e from 0.08 to 0.062. The parameter stays in the range measured experimentally [30,31]. Interestingly, the more gradual dependence is also used in dynamic clamp experiments [9,13]. As shown below, the modification is very important because it is responsible for the high frequency firing achieved during somatic NMDA receptor activation. The calcium equation (2) represents balance between Ca 2+ entry via the L current and Ca 2+ removal via a pump. After entering the cell, Ca 2+ binds to a buffering protein. Buffering is assumed instantaneous and is taken into account by multiplying by the buffering coefficient b. This coefficient is the ratio of free to total intracellular Ca 2+ . Therefore, [Ca 2+ ] represents free intracellular Ca 2+ concentration, and the parameters were estimated before in experiments [20].

Morphologies
The biophysical mechanisms (currents and pumps) described above were inserted into two different morphologies: a single compartment, and a reconstruction of a DA neuron [20]. The equations introduced above fully describe the one-compartmental model. Computer simulations of this model were done using XPPAUT [32].
Simulations in the reconstructed morphology [20] were done using NEURON (http://www.neuron.yale.edu). All biophysical mechanisms were distributed homogeneously along the dendrites and in the soma. The dendrites and soma were divided into 50 mm-long equipotential segments. The value of 60 ohm|cm was used for the cytoplasmic resistivity, which is in the range of previously used (see e.g. [33]).
The models are published online in Model-DB.

Reconstructed Morphology Model
Somatic NMDAR activation evokes high-frequency firing and subthreshold oscillations. Experiments [9,13] show that NMDA receptor activation restricted to the soma effectively evokes high-frequency oscillations. However, modeling reproduced frequency elevation under dendritic, but not somatic, NMDA synaptic stimulation [12]. We correct this discrepancy: Fig. 2 A shows the effect of simulated somatic NMDA receptor activation implemented in a reconstructed morphology of the DA neuron. The NMDAR current is injected into the soma at 1000 ms for the duration of 1000 ms. The activation and removal of the current are modeled as instantaneous step functions and cause rapid transition between high-and low-frequency oscillations. The frequency reaches 25 Hz. A similar frequency increase occurs if a dendritic, but not somatic, or both the dendritic and somatic NMDA receptors are activated (Fig. 2B, D). An excessive NMDAR conductance results in depolarization block (Fig. 2 C).
For the whole-neuron NMDAR activation (bath NMDA), the frequency above 20 Hz (typical for in vivo bursts) is achieved at much lower NMDAR conductance levels ( g g NMDA &7mS=cm 2 ). These levels are well within the physiological range. Several times higher NMDAR conductance is necessary to achieve that frequency for the soma-only excitation ( g g NMDA w22mS=cm 2 ). This is consistent with experiments [13]. The frequency rise depends on the spread of NMDA receptor activation as well as the focus of its application. We simulated iontophoresis of an NMDAR agonist to the distal dendrites of the neuron. No high frequency was achieved when any single branch of the dendritic tree was stimulated. The oscillations took a complex form of intermixed singular high-and multiple lowamplitude oscillations (mixed mode; data not shown). The total surface area of the stimulated dendrite was similar to, or even greater than the somatic one. Only when three or more distal dendrites were stimulated together (Fig. 3A, red), did the frequency rise significantly. Compared to the soma, approximately three times greater area of distal dendrites had to be affected by NMDA to evoke the high frequency. This is because distal dendrites need to synchronize across the dendritic tree for their concerted influence to dominate the soma and proximal dendrites (Fig. 3B). When there are stimuli that partition the dendritic tree, it is no longer synchronized. In Fig. 3B, even the dendrites that receive the NMDAR stimulation are not synchronized one-to-one: the dendrite 1 files an extra spike per the period (red trace). This spike is not enough to evoke a spike in the soma, and only when all three dendrites fire, the somatic spike follows. No further increase in the NMDAR conductance synchronizes the dendrites one-toone. As a result, distal dendritic NMDAR stimulation is less effective than focused proximal or somatic.
In the DA neuron, the oscillations persist in the subthreshold range after the blockade of the spike-producing currents (TTX application). Many publications document this for the lowfrequency background activity (see e.g. [20]), and this was recently extended to the NMDA-evoked high-frequency activity [9]. Fig. 2D also shows the frequency responses during simulated blockade of the spike-producing fast sodium current (thin curves). The amplitude of oscillations becomes much smaller (see Figs. [4][5][6][7][8][9], whereas the frequency does not change significantly in the model. These matching frequencies in the model with and without the spiking currents mirror those in experimental data for TTX vs. control [9,20]. The model does not allow for the comparison of the regularity of spiking and subthreshold oscillations [17] because no noise is included and the oscillations are perfectly regular. The match between spiking and subthreshold oscillations justifies our first reduction: in the rest of the manuscript, we show results for the model with no sodium spike-producing currents. The nonspiking model has the minimal number of biophysical components and the minimal number of variables (three: v, [Ca 2+ ] and n) that are required for the described phenomena.
Tonic AMPAR stimulation and applied depolarization lead to early blockade of oscillations. By contrast to NMDA, applied depolarization or tonic AMPAR synaptic stimulation is unable to elevate the frequency in the DA neuron to the levels observed during bursting [9,10,11]. In our model, tonic applied depolarization elevates the frequency to 6 Hz only and then leads to blockade of oscillations ( Fig. 4 A, C). Tonic activation of the AMPAR current also led to blockade of oscillations rather than a frequency increase (Fig. 4 B, D). The range of the AMPAR conductance where oscillations persist is strikingly smaller than that for the NMDAR conductance. Note that complex, mixedmode-type oscillations are displayed during and after AMPAR stimulation. This is very similar to voltage traces obtained after the stimulation in experiments ( [9], Fig. 4A) Therefore, the frequency remains limited below 10 Hz under both AMPAR activation and applied depolarization, whereas NMDAR activation breaks this frequency limitation.
A combination of calcium-and voltage-dependent currents repolarizes the cell. Two lines of evidence suggest that the fast currents hyperpolarizing after a spike (fast AHP complex), including the delayed rectifier potassium current, are weak in the DA neuron. First, our modeling shows that applied depolarization would evoke much higher frequencies if the delayed rectifier was stronger. This would erase the difference between the responses to AMPAR and NMDAR synaptic stimulations. Second, spiking becomes unreliable and the neuron becomes prone to the cessation of oscillations when the SK-type calciumdependent potassium current is blocked pharmacologically. The blockade was used in vitro to elicit bursting [34], but within or without the bursts, spiking quickly collapsed to a plateau.
On the other hand, the neuron demonstrates calcium-independent firing, as we mentioned in the Introduction. This requires a calcium-independent repolarizing mechanism. We have introduced a voltage-dependent potassium current to play this role. The calibration of the current resulted in the half-activation around 250 mV, slow activation (,62 ms) and even slower deactivation (,362 ms). The slowness of the current, and especially that its deactivation is much slower than activation, suggests the ERG current as the best candidate [35]. The current sustains pacemaking in the absence of the SK current, but only in a very small interval of applied depolarization (Fig. 5 A&B). We found that this property is achieved by a close match between the half-activations of the Ca 2+ current and our ERG current in the model. Our nullcline analysis below explains this connection. The slowness of the current determines that only moderate increase in the frequency of background firing is observed upon the blockade of the SK current (Fig. 5A), which matches experimental data [36,37]. The difference between the activation and deactivation timescales allows the ERG current to maintain the voltage waveform with a narrow peak and slow release from the hyperpolarized state. Additionally, we slow down the activation of the current to make it negligible in the presence of the SK current (Fig. 6). This makes the SK current dominant in the oscillatory cycle. Upon the blockade of the SK current, the switch from it to the ERG current determines the narrowing range of applied depolarization where the neuron sustains oscillations.
Our simulations predict that the NMDA-evoked high-frequency oscillations are robust when the SK current is blocked. This is in contrast to the property that oscillations become even more fragile in response to applied depolarization (Fig. 5F). Oscillations persist in a wide range of NMDAR conductance density and display elevated frequencies (Fig. 5C,E). Moreover, moderate NMDAR stimulation allows oscillations to persist at much higher applied depolarization and rescue the neuron form the cessation of oscillations (Fig. 5D). This prediction extends the observation that NMDA-evoked high-frequency firing persists after blockade of the SK current by apamin [9]. In Supplement S1, we illustrate both low and high-frequency firing during the blockade of the SK current in the full model with spike-producing currents (Fig. S1B in Supplement S1). Thus the ERG current supports and amplifies the distinction between the responses to applied depolarization and NMDAR stimulation.

One-compartmental Model
Minimal representation of the DA neuron. Next we show that a one-compartment representation is sufficient for reproducing distinct frequencies of subthreshold oscillations evoked in response to NMDA vs. AMPA receptor stimulation and applied depolarization. The above model implements minimal biophysical properties, but the number of compartments representing the reconstructed morphology makes it very complex. A reduction is essential for explaining our results. If all parts of the neuron are in synchrony, a one-compartment model will work. Fig. 7 shows typical responses of the one-compartment model to AMPAR and NMDAR activation and somatic depolarization. The results of stimulations closely match those for the reconstructed morphology and justify the assumption of synchrony (we compare the single compartment to the bath application of the agonists above). Now the radius of the compartment in Eq. (2) does not correspond to any particular part of the neuron any more. It represents an ''average'' compartment because all properties, including the frequency are effectively averaged by the coupling when the compartments are in synchrony. By choosing a value for the radius, we calibrated the model to reproduce maximal frequencies similar to those obtained during applied depolarization and NMDAR activation in experiments. Interestingly, the radius takes the value of 0.65 mm, which corresponds to fine distal dendrites. To understand this, let's think of a reduced model with the same number of compartments, but where the radii of compartments are all the same and have the diameter of a distal dendrite. Thus, we replaced all proximal dendrites and the soma with distal dendrites. We know that such compartments synchronize as one and produce a mode similar to the reconstructed morphology. This suggests that the total contribution of the distal dendrites to the overall transmembrane current is very substantial, and distal dendrites control the frequency of the whole neuron. Only when there is a stimulus that splits the dendritic tree, the soma and proximal dendrites receive the opportunity to contribute more. The single-compartment model does not account for such cases, but it reproduces all other properties described above.
Mechanism for the high frequency: flattening the voltage nullcline. To explain the mechanism underlying the frequency responses, we simplify the model by blocking our ERG current. As a result, the model is reduced to a system of two variables: the voltage and calcium concentration (Eq. 1&2). Fig. 8A shows time series for these two variables at the onset of NMDA-evoked highfrequency oscillation. To show a good correspondence with the reconstructed morphology model that includes the spike-producing currents, we show the same transition in Fig. S1A of Supplement S1. In Fig. 8A, at the onset, the amplitude of Ca 2+ oscillations is dramatically reduced. This underlies the frequency increase. A less steep voltage dependence of the NMDA current compared to the Ca 2+ current determines the decrease in the amplitude of Ca 2+ concentration. To explain this, we introduce a simple mathematical tool.  A basic mathematical formalism allows for geometric visualization of the changes introduced by NMDAR activation. Plotting the voltage and Ca 2+ concentration against each other at every instance of time results in a cycle that represent oscillations (Fig. 8B). We investigate the dynamics by plotting nullclines of the system. These are two curves on which either Ca 2+ or voltage equilibrates ( d½Ca 2z dt = 0 or dv dt = 0, respectively). Oscillations circumscribe folding of the voltage nullcline as shown in Fig. 8B. The intersection of the two nullclines is an equilibrium state of the model. This state corresponds to the blockade of oscillations. It always exists, but may be unstable (repulsive) and, consequently, not displayed by the neuron. The oscillations are impossible if the equilibrium state interrupts the oscillatory cycle (Fig. 8C). The voltage dependence of the NMDAR current is central to shaping the voltage nullcline. The changes to the nullclines under increasing maximal NMDAR current density and applied depolarization are shown in Fig. 9A. During NMDAR activation, the intersection of the nullclines remains between the extrema and remains an unstable equilibrium. During a current injection, the  minimum of the voltage nullcline shifts across the calcium nullcline. As a result, the equilibrium state at their intersection interrupts the oscillatory cycle, and oscillations die out. The equilibrium becomes stable (attractive). This explains why, oscillations die out with a small elevation in the applied current, but not at the higher NMDAR current density.
Additionally, the NMDAR activation makes the voltage nullcline more flat. More precisely, the extrema become significantly closer to each other along the Ca 2+ concentration axis, whereas they stay well separated along the voltage axis (Fig. 9A). Oscillations span the interval of the Ca 2+ concentration between the extrema. The length of this interval is the primary determinant of the period of oscillations because Ca 2+ concentration changes slowly. NMDAR shortens the interval of Ca 2+ concentration between the extrema. As a result, NMDAR significantly reduces the period by reducing the amplitude of Ca 2+ oscillations. Taken together, the high-frequency oscillations are achieved by flattening, but not unfolding the voltage nullcline.
Extracellular magnesium blocks the NMDA receptor at a low voltage. This defines the NMDAR current half-activation, and the higher the magnesium concentration, the higher the half activation. The dependence of the frequency on Mg 2+ concentration was measured in dynamic clamp experiments [9]. Our model reproduces this dependence (Fig. 10). At a higher Mg 2+ concentration, NMDAR current remains constantly blocked, and the frequency remains low. At a very low Mg 2+ concentration, the NMDAR current is always open and looses its voltage dependence. It acts as the AMPA receptor and suppresses oscillations. In the middle range, the frequency increases with decreasing Mg 2+ concentration and reaches above 100 Hz in its maximum. The properties closely match the experimentally obtained dependence and expand it to the higher Mg 2+ concentration. Note that in [9] a different coefficient was used in the formula for the NMDAR current dependence, and our Mg 2+ concentration must be divided approximately by the factor of 3 to obtain the same scale.
The other parameter of the dependence, m e defines its slope. This slope determines how the voltage nullcline transforms as the NMDAR current is introduced. The slope of the voltage nullcline is directly related to the slope of the I-V curve for the neuron. In particular, a negative resistance caused by regenerative opening of depolarizing channels corresponds to the middle part of the voltage nullcline, where it has a positive slope. This part is steeper if the activation curves for the currents are steeper. Conversely, flattening of the voltage nullcline (Fig. 9A) reflects a gradual voltage dependence of the NMDAR current (Fig. 1B). The increase in the frequency corresponds to the gradual flattening of the voltage nullcline as the NMDAR contributes a greater portion of the depolarizing current. For comparison, Fig. 9B shows the nullclines for the NMDA current with a steeper voltage dependence (m e~0 :08). The vertical distance between the extrema of the voltage nullcline remains large. Accordingly, the frequency of oscillations remains low for any conductance of the NMDAR current. Therefore, the ability of the NMDAR current to evoke the high frequencies is linked to the moderate slope of its voltage dependence. This gives an easily testable prediction: a virtual NMDA current with a steeper voltage dependence will not elevate the frequency.
Dynamical consequences of combining the repolarizing currents. The inclusion of the two repolarizing currents makes the model three-dimensional, where each of them is capable of sustaining oscillations. The currents form two negative feedback loops. They are shown in Fig. 11, where a hammerhead means inhibition and an arrow means activation. However, it's impossible to simultaneously observe two different rhythms generated by the separate feedback loops. This is because the loops share a component, which is the voltage. Any manipulation that breaks a connection breaks at least one feedback loop. This distinguishes the structure of the model from coupled oscillators. Such structure is well known for the circadian clock and was called interlocked feedback loops.
Voltage oscillations of a complex form are often observed in experiments (see e.g. [9]), and our model suggests a way they can be generated. The one-compartment model displays bistability, in which stable oscillations coexist with a stable equilibrium state. In this case, the observation of oscillations or their blockade is dependent on the initial conditions (e.g. holding potential). Subtle differences in experiments may modulate initial conditions and lead to different coexisting modes. In Fig. 12A, two types of oscillations very different by their amplitude are observed in the model in the absence of stimulation. The high-amplitude oscillatory solution corresponds to the time series in Fig. 7 without stimulation (the first and the last second of each time-series). None of these stimulations switches the model into the low-amplitude mode shown in Fig. 12. Only after a stronger applied current or AMPAR stimulation that suppresses oscillations, the one-compartment model remains close to the equilibrium state and displays the low-amplitude oscillations. NMDAR stimulation also removes bistability, but switches the model into the high-amplitude mode. Interestingly, even after the suppression of oscillations achieved by excessive NMDAR stimulation, as soon as the stimulation ends, the model returns to the high-amplitude mode. The latter holds in the reconstructed morphology model (Fig. 2 C), but the model returns to the high-amplitude mode after suppression of oscillations by AMPAR stimulation or applied current, in contrast to the single compartment. Our simulations show different mixed modes and long complex transients after the stimulation (see Figs. 2 B&C, 4 A&B). We may view the mixed modes as switching between the Figure 10. The dependence of the oscillation frequency on the magnesium concentration at a fixed NMDAR current density ( g g NMDA~1 0mS=cm 2 ). doi:10.1371/journal.pone.0069984.g010 high-and low-amplitude modes by interaction among compartments. Thus, we suggest that complex oscillatory modes in the reconstructed morphology are based on bistability in the minimal single-compartment model. Multiple mixed mode solutions coexist with each other. Thus, multystability and complex modes displayed in our model may explain apparently discrepant experimental results obtained for the DA neuron in different groups (e.g. [9], [11]).

The Role of the DA Neuron Morphology
We have made two conclusions on the role of the DA neuron morphology: First, the morphology does not determine the capacity of the DA neuron to generate a high-frequency firing differentially in response to NMDAR, but not AMPAR stimulation or applied depolarization. Second, both the spread and the focus of NMDAR stimulation along dendrites affect the transition to the high-frequency firing. Both statements and a relation between them deserve attention.
Puzzling firing properties of the DA neuron trigger a search for very complex mechanisms, including its morphology and heterogeneity of membrane currents. In our previous study [12], the interaction of the dendrites and the soma determines the firing rate. However, this mechanism was disproved in experiments [9,13]. The new property completely incompatible with the old mechanism was the ability of the focal somatic NMDA stimulation to elicit the high-frequency firing. This is the first property we have reproduced in our new model. The model is presented in two morphologies: (1) a reconstruction of a DA neuron, and (2) a single compartment that ignores the spatial structure of the neuron. A great correspondence between the two representations suggests that an unusual morphology of the DA neuron [38] is not necessary for the neuron to distinguish between NMDA and AMPA receptor stimulation. For this purpose, the whole dendritic tree may act as one compartment. On the other hand, the dendritic morphology plays a role in synaptic integration of a heterogeneous stimulation. By the comparison of the two representations of our model, we made the conclusion that distal dendrites control the frequency of the neuron, but only when the dendrites are synchronized. A focal dendritic stimulation separates a part of the dendritic tree and impedes synchronization. In this case, the single compartment representation is not correct. The morphology determines the minimal portion of the dendritic tree that must receive NMDA stimulation to evoke high frequency firing. That was approximately three times greater area of distal dendrites compared to the soma in our model. By contrast, the somatic NMDAR stimulation spreads to all dendrites evenly and does not violate their synchronization. This leads to a surprising conclusion: the soma is at a much better location to manipulate the firing frequency, even though the frequency is controlled by distal dendrites.

Multiple Pacemaking Mechanisms
A number of recent publications challenge the calciumpotassium pacemaking mechanism in the DA neuron and suggest Na + -based pacemaking instead [11,[16][17][18]. The Na + -based pacemaking requires a calcium-independent repolarizing current. Several fast voltage-gated hyperpolarizing currents have been identified in the neuron [39]. However, the fast AHP currents must be weak in he DA neuron because otherwise it would make the responses to AMPA and NMDA receptor activation indistinguishable and firing in the presence of apamin perfectly robust. Therefore, a slow voltage-sensitive hyperpolarizing current is necessary to sustain calcium-independent pacemaking, and no study suggests what current may play this role. We have introduced a simple voltage-dependent current and its calibration suggests the ERG current for this role. ERG1 proteins are expressed in the substantia nigra pars compacta neurons [23]. A slow AHP similar to that produced by the ERG current was measured in the DA neuron [24], although the current has not been fully characterized. Early studies of the ERG current showed its deactivation to be much slower than activation [35]. This is exactly what our calibration gives: the activation in the range of tens of msec. and deactivation in the range of hundreds of msec. Our half-activation of 247.4 mV closely match the experimental estimates. The slope of the activation function measured in Purkinje neurons is 5 mV, and we use 2 mV in our model. The difference may be explained by the variability of the current in different neuron types [35,40]. The ERG current is suggested to be involved in termination of bursts [24,41] and affect pacemaking in the DA neuron [25]. Our study attributes a new major role to this current for repolarizing in calcium-independent firing.
Oscillatory mechanisms in the DA cell have been intensely discussed. Originally, Ca 2+ and Ca 2+ -dependent potassium currents were shown to sustain subthreshold oscillations as well as spiking in the DA neuron (see e.g. [19][20][21]). Then, a significant contribution of Na + currents to pacemaking was shown [11,[16][17][18]42], and the mechanism that generates spiking was proposed to be different from one generating the subthreshold oscillations [17]. Then the conclusion that the Ca 2+ and Na + currents contribute to the same oscillatory mechanism has been reached [22]. Thus, spike-producing currents alone are unable to sustain firing and require subthreshold drive, consistent with our conclusion above. However, this mechanism cannot sustain Ca 2+ -independent pacemaking found in the DA neuron [11,[16][17][18]. We propose the second subthreshold oscillatory mechanism that can play this role. The two mechanisms are defined by two negative feedback loops (Fig. 11). A Ca 2+ -dependent and a voltage-dependent repolarizing current constitute the loops, and each of the loops can sustain oscillations. The voltage-dependent loop represents the Ca 2+ -independent pacemaking. In our model, depolarization is still provided by the Ca 2+ current, but once the SK current is blocked, the Ca 2+ specificity of the current is irrelevant because Ca 2+ concentration does not affect the voltage. In accordance with that, the voltage dependence, but not Ca 2+ specificity of the Ca 2+ current has been found critical for pacemaking [13]. Thus, our model supports the idea that there are multiple pacemaking mechanisms in the DA neuron. Our prediction is that the mechanisms must sustain subthreshold oscillations, which then drive spiking.
Another prediction can be tested in dynamic clamp experiments. Not only blockade of the ERG current, but also its enhancement does not significantly change the background firing in the presence of the SK current. When the SK current is blocked, our model predicts that a small enhancement of the ERG current blocks the oscillations. By contrast, a significant reduction of the ERG current (4.8 to 0.9 uS/cm 2) only changes the waveform: the voltage stays longer in the depolarized phase, which was observed in experiments and called plateau potentials [41]. Another predictions can be tested by a simple current clamp experiment when the ERG current is blocked and the SK current stays intact. Our simulations show that oscillations persist at a much higher applied depolarization (data not shown) in this case than when the ERG current is present. Our explanation is that, during applied depolarization, the ERG current develops a constant component that shunts the active currents. The mechanism is similar to shunting inhibition [43].
The structure of the model (Fig. 11) links firing patterns of the DA neuron with research in another branch of biology and may explain the role of this structure. Such structure defines oscillatory mechanisms that underlie the circadian rhythm [44]. This structure is called interlocked feedback loops. The loops are interlocked in the sense that it is impossible to separate them without losing the oscillations. This interlocked feedback loop structure has been suggested to improve robustness of the circadian rhythm [44]. Accordingly, blockade of the SK current in the DA neuron, which brakes one of the feedback loops, dramatically reduces robustness and regularity of firing [36,37]. Interaction of multiple ion channels has been suggested to ensure robustness of pacemaking with respect to perturbations in individual components [17]. Robustness of oscillations can be studied using minimal models, in which biophysical details are omitted and only the structural properties are reflected. These models are the same for different applications, such as regulatory networks or electrical activity of a neuron. Therefore, the interlocked feedback loop representation provides an interdisciplinary framework for the analysis of robustness.

The Role of NMDA in Eliciting Bursting
A number of studies associate bursting specifically with NMDA receptor activation both in vivo and in vitro [5,7,9,45,46]. The distinction between AMPA and NMDA receptor currents that allows the DA neuron to respond differentially is the NMDAR current voltage dependence [9,13]. Our study supports the importance of the voltage dependence of the NMDAR current. We have shown that NMDAR and AMPAR activation elicit distinct responses because they interact differently with the mechanism for pacemaking in DA neurons. The NMDAR current expands the interval of negative resistance by its regenerative opening. As a result, NMDAR effectively counteracts decreasing amplitude of oscillations with respect to the voltage. Simultaneously, it significantly reduces the amplitude of the Ca 2+ oscillations. This leads to a pronounced elevation of the frequency because Ca 2+ concentration is a slow variable and its variations determine the period. This requires the slope of the voltage dependence of the NMDAR current to be moderate. We have shown that NDMAR current with a steeper voltage dependence never elevates the frequency. On the other hand, the NMDAR current with a gradual voltage dependence cannot restore regular pacing when the Ca 2+ current is blocked in the model, which replicates experiments with nimodipine-silenced DA neurons [13]. The ability of the Ca 2+ current, but not NMDA receptor current to sustain pacemaking was attributed to a difference in their halfactivation. We stress that the slope of the voltage dependence plays even more important role: the frequency elevation would not be possible with a steeper dependence at any half-activation. Therefore, the model predicts that the ability of the NDMAR current to evoke the high frequencies is determined by the moderate slope of its voltage dependence.
Our results do not exclude the involvement of AMPA receptors with bursting in DA neurons. Previous experimental [11] and modeling [47] studies show that both AMPA and NMDA receptors contribute to burst firing. Here, the major distinction with our work is that the receptors were activated with highfrequency trains of electrical stimulation. Such stimulation evokes short powerful AMPA EPSPs by synchronous activation of the receptors. Our model is capable of following a pulsatile AMPA stimulation, but shows that stronger bursting results from an NMDA stimulation of the same intensity (Fig. S2 of Supplement S1). Similarly, pulses of applied depolarizing current at 25 Hz separated by weak hyperpolarizations evoke spiking at this frequency [9]. Our model replicates this experimental finding (data not shown). An asynchronous AMPA tone that sums into a plateau potential fails to elevate the frequency [9]. This suggests that in the in vivo conditions of the previously published studies [5][6][7] AMPA receptors activate asynchronously.
Our simulations show that NMADR stimulation dramatically increases the basal level of Ca 2+ concentration (see Fig. 8A). This increase is due to the opening of the calcium current that follows depolarization. The flux of calcium via L-type Ca 2+ channel is proposed to be responsible for the oxidative stress in the DA neuron and cause degeneration of the neurons leading to Parkinson disease [16]. Calcium influx through NMDA receptors will also contribute to the neurotoxicity [48]. Thus, our study suggests a pivotal role of NMDAR activation in degeneration of DA neurons.

Connections to other Modeling Studies
The electrophysiology of the DA neuron has been modeled by several groups, and each model was focused on a particular property. The calcium-potassium mechanism for background firing was first modeled in a single compartment by [49], and then extended throughout the length of the dendrites by [20]. This gave rise to a coupled oscillator representation of the DA neuron. Here, we start with this representation and show that it can be reduced back to a single compartment.
Rhythmic bursting evoked by bath application of NMDA agonist in DA neurons was modeled in a number of studies [29,50,51]. These studies were focused on a specific feature for this type of bursting -the rhythmic plateau potentials that underlie bursts. As we emphasize above, the high-frequency firing cannot be a simple consequence of a plateau potential in DA neurons, and the present study offers a unique mechanism for the intraburst firing, which can be combined with plateau potential generation.
Bursting achieved by the blockade of the SK type Ca 2+dependent potassium current was modeled in a number of papers [22,41,47,49,52]. The ERG current was suggested to provide extraburst hyperpolarization in the absence of the SK current [25,41]. We show that, surprisingly, such a slow current can sustain fast firing within the bursts. A stronger activation of this current can also end the burst and provide extraburst pauses, as shown earlier. The intraburst high-frequency firing was assumed to be a simple consequence of the plateau potentials in most of these papers, except for two. [52] show clustering spikes into highfrequency bursts only in the presence of an SK current blocker. NMDA stimulation alone could not evoke bursting in their model, and our model covers this gap. The authors of [22] mostly focus on the combination of Na + and Ca 2+ depolarizing currents and managed to reconcile apparently discrepant experimental results on their contribution (see e.g. [17]). However, their model continues firing in the absence of the SK current due to additional Ca 2+ -dependent repolarization by the Ca 2+ pump, whereas our model offers a Ca 2+ -independent repolarizing current, which works for Na + -based pacemaling. Both cases are possible, and which repolarizing mechanism works when the SK current is blocked is a question for future experiments. While [22] suggests that the amplitude of calcium oscillations increases upon blockade of the SK current, our study predicts that the whole range of calcium concentrations shifts up more significantly. Thus, using this prediction, a future Ca 2+ -imaging experiment will determine the mechanism of repolarization.
Finally, the mechanism of high-frequency oscillations presented here is similar to one presented in our recent publication [53,54]. Here and in these publications, flattening of the voltage nullcline determines the frequency elevation. The major difference is that here we link this flattening to the gradual dependence of the NMDA current. The alternative way presented before was based on the stretching of the voltage nullcline along Ca 2+ concentration by the nonlinear dependence of the SK current. This stretched shape was initially formed so that a stimulation could shift the range of Ca 2+ concentration to a linear part of the SK current dependence, and the nullcline became compressed. The elevation of frequency obtained this way was more steep compared to one obtained here (Fig. 2). Interestingly, the function of the NMDAR current voltage dependence used in dynamic clamp experiments had the reduced slope that we use in this work [9,13]. This favors the mechanism presented here, but does not exclude the possibility that both mechanisms are combined in the neuron.

Author Contributions
Conceived and designed the experiments: AK JH. Performed the experiments: AK JH. Analyzed the data: AK JH. Contributed reagents/ materials/analysis tools: AK JH. Wrote the paper: AK JH.