Calcium Dynamics Underlying Pacemaker-Like and Burst Firing Oscillations in Midbrain Dopaminergic Neurons: A Computational Study

B. Amini,1 J. W. Clark, Jr.,1 and C. C. Canavier2

 1Department of Electrical and Computer Engineering, Rice University, Houston, Texas 77005; and  2Department of Psychology, University of New Orleans, New Orleans, Louisiana 70148


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Amini, B., J. W. Clark Jr., and C. C. Canavier. Calcium Dynamics Underlying Pacemaker-Like and Burst Firing Oscillations in Midbrain Dopaminergic Neurons: A Computational Study. J. Neurophysiol. 82: 2249-2261, 1999. A mathematical model of midbrain dopamine neurons has been developed to understand the mechanisms underlying two types of calcium-dependent firing patterns that these cells exhibit in vitro. The first is the regular, pacemaker-like firing exhibited in a slice preparation, and the second is a burst firing pattern sometimes exhibited in the presence of apamin. Because both types of oscillations are blocked by nifedipine, we have focused on the slow calcium dynamics underlying these firing modes. The underlying oscillations in membrane potential are best observed when action potentials are blocked by the application of TTX. This converts the regular single-spike firing mode to a slow oscillatory potential (SOP) and apamin-induced bursting to a slow square-wave oscillation. We hypothesize that the SOP results from the interplay between the L-type calcium current (ICa,L) and the apamin-sensitive calcium-activated potassium current (IK,Ca,SK). We further hypothesize that the square-wave oscillation results from the alternating voltage activation and calcium inactivation of ICa,L. Our model consists of two components: a Hodgkin-Huxley-type membrane model and a fluid compartment model. A material balance on Ca2+ is provided in the cytosolic fluid compartment, whereas calcium concentration is considered constant in the extracellular compartment. Model parameters were determined using both voltage-clamp and calcium-imaging data from the literature. In addition to modeling the SOP and square-wave oscillations in dopaminergic neurons, the model provides reasonable mimicry of the experimentally observed response of SOPs to TEA application and elongation of the plateau duration of the square-wave oscillations in response to calcium chelation.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Midbrain dopaminergic neurons (DA) have been implicated in several important functions in humans, including movement, attention, learning, and reinforcement as well as in the etiology of Parkinson's disease and various thought disorders, such as schizophrenia, attention deficit disorders, and depression (Carlson 1992). Therefore it is essential to understand how the firing pattern in these neurons is generated and modulated. Whereas in vivo DA neurons can exhibit either single-spike or burst-firing activity as well as spontaneous shifts between modes (Freeman et al. 1985; Grace and Bunney 1984), in a slice preparation a regular, pacemaker-like firing pattern usually is observed, presumably due to the loss of synaptic afferents. In this study, we will focus on the intrinsic oscillatory activity of DA neurons in a slice preparation.

In addition to the endogenous regular single-spike firing exhibited in a slice preparation, burst firing has been induced by the application of apamin (Ping and Shepard 1996) or N-methyl-D-aspartic acid (NMDA) (Johnson et al. 1992). The firing pattern in these neurons is particularly significant because burst-firing in vivo is associated with behaviorally relevant appetitive and sometimes novel stimuli and possibly with movement (for a review, see Overton and Clark 1997). This study will address the biophysical basis for the calcium-dependent oscillations underlying regular single-spike firing and apamin-induced burst firing.

Because midbrain DA neurons have an unusually depolarized firing threshold (-33 ± 15 mV) (Grace and Onn 1989), the membrane must be depolarized repeatedly from rest (-60 mV) to this threshold to fire repetitively. Indeed, such an underlying slow oscillatory potential (SOP) has been observed when the sodium-mediated spikes are blocked by tetrodotoxin (TTX) and is enhanced by the application of the potassium channel blocker tetraethyl ammonium (TEA). The depolarizing phase of the SOP has been variously termed the pacemaker-like slow depolarization (PLSD) or simply the pacemaker-like depolarization (PLD), which has been shown to be calcium-dependent.

The SOP sometimes can be converted into a square-wave oscillation by the application of apamin (Nedergaard et al. 1993; Ping and Shepard 1996), revealing a crucial role for the apamin-sensitive current IK,Ca,SK in the repolarization of the SOP. Just as the SOP is the underlying oscillation that sets the rhythm for pacemaker-like single-spike firing, the square-wave oscillation sets the timing for burst firing; in the absence of TTX, a burst of spikes is generated on the plateau of the square-wave and quiescence during the trough (Gu et al. 1992; Ping and Shepard 1996; Shepard 1993; Shepard and Bunney 1988). In some cases, apamin induced irregular firing instead of burst firing (Gu et al. 1992; Ping and Shepard 1996), and in those cases, the subsequent application of TTX and TEA revealed irregular high-threshold calcium spiking (Ping and Shepard 1996). The generation of irregular firing patterns, however, is not addressed in this study. Neither the SOP nor the square-wave should be confused with the slow oscillation in membrane potential revealed by TTX in the presence of NMDA, which also can underlie burst firing, appears to be sodium rather than calcium dependent (Johnson et al. 1992) and is also not addressed in this study.

The goal of this study was to test hypothesized mechanisms for the SOP and the square wave. On the basis of data that suggest that dihydropyridines abolish the SOP (see DISCUSSION), we hypothesize that the SOP is generated by interplay between the L-type calcium current ICa,L and the calcium-activated IK,Ca,SK. We further hypothesize that in the presence of apamin, the calcium-mediated inactivation of ICa,L is responsible for repolarization of the square wave. We also hypothesize that the mechanisms driving these phenomena are located in or near the soma. Hence we constructed a quantitative somatic model based on voltage-clamp, morphological, and calcium-imaging data to confirm that the currents described under voltage clamp are capable of generating the electrical activity observed under current clamp. Minimum expectations for model performance included production of SOPs similar to those of Ping and Shepard (1996) in the presence of TTX and TEA, SOPs under conditions where either the calcium current ICa,T or ICa,N were blocked but not when ICa,L is blocked, and square waves similar to those observed (Ping and Shepard 1996) in the presence of apamin.


    MODEL DEVELOPMENT
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Our model of the DA neuron consists of a single compartment Hodgkin-Huxley (HH)-type parallel conductance membrane model (Fig. 1A) and a lumped fluid compartmental model (Fig. 1B). The HH equivalent circuit is composed of a somatic membrane capacitance of 15.8 pF (Kang and Kitai 1993b), shunted by resistive ion-selective channels, as well as ion pumps, a sodium-calcium exchanger, and other ionic currents. The lumped fluid compartment model (Fig. 1B) consists of an intracellular compartment containing constant concentrations of Na+ and K+, and a calcium buffer (presumably calmodulin). A material balance for Ca2+ describes the time rate of change in cytosolic Ca2+ concentration. The extracellular space is assumed to have a relatively large volume, so that the ionic concentrations of Ca2+, Na+, and K+ there are assumed to be constant.



View larger version (21K):
[in this window]
[in a new window]
 
Fig. 1. Diagram of the dopaminergic (DA) neuron cell model. A: electrical equivalent circuit for the membrane; ICa is the combined calcium current consisting of ICa,T, ICa,N, ICa,L, and ICa,HVA. IK,tot is the combined voltage-dependent potassium current consisting of IK and IA.B: fluid compartmental model, including intracellular and extracellular spaces.


    MEMBRANE CURRENTS
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Under space-clamp conditions, the differential equation describing the time-dependent changes in the membrane potential (V) is
<IT><A><AC>V</AC><AC>˙</AC></A></IT><IT>=</IT>−(<IT>I</IT><SUB><IT>Ca,T</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>Ca,L</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>Ca,N</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>Ca,HVA</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>K</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>A</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>h</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>Na,K</IT></SUB> (1)

<IT>+</IT><IT>I</IT><SUB><IT>Ca,pump</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>Na,Ca</IT></SUB>)<IT>/</IT><IT>C</IT><SUB><IT>m</IT></SUB>
where Cm is the whole cell membrane capacitance. Model currents include a T-type calcium current (ICa,T), an L-type calcium current (ICa,L), an N-type calcium current (ICa,N), a residual high-voltage activated calcium current (ICa,HVA), a delayed rectifier current (IK), a transient outward current (IA), a small conductance calcium-dependent potassium current (IK,Ca,SK), a hyperpolarization activated current (Ih), a sodium-potassium pump current (INa,K), a calcium pump current (ICa,pump), and a sodium-calcium exchanger (INa,Ca).

In the current descriptions that follow, the HH-type activation and inactivation gating variables are solutions of the familiar first-order differential equations described as
<IT><A><AC>z</AC><AC>˙</AC></A></IT>(<IT>V</IT><IT>, </IT><IT>t</IT>)<IT>=</IT><FR><NU><IT><A><AC>z</AC><AC>&cjs1171;</AC></A></IT>(<IT>V</IT>)<IT>−</IT><IT>z</IT>(<IT>V</IT><IT>, </IT><IT>t</IT>)</NU><DE><IT>&tgr;</IT><SUB><IT>z</IT></SUB>(<IT>V</IT>)</DE></FR> (2)
where <A><AC>z</AC><AC>&cjs1171;</AC></A>(V) is the steady-state value of the general gating variable z at membrane voltage V. We have characterized the steady-state gating variable <A><AC>z</AC><AC>&cjs1171;</AC></A>(V) by a sigmoidal, or Boltzman-type relationship, and the time constant tau z(V) by a Gaussian relationship.

Individual ionic membrane currents were characterized by fits to published voltage clamp experiments where available. The temperature at which the experiments were performed ranged between 30 and 35°C, thus temperature adjustments (e.g., Q10) were not needed during final integration of currents into the model. To obtain fits to the currents, membrane potential was held constant and the differential equations characterizing the gating variables for each current were integrated numerically and substituted into the appropriate equation for ionic current. Examples of fits to the component currents of the model are shown in Fig. 2.



View larger version (32K):
[in this window]
[in a new window]
 
Fig. 2. Model-generated approximations to voltage-clamp data for several membrane currents. Steady-state activation and inactivation characteristics of the currents were characterized by fits (---) to experimental data (- - -). A: ICa,T, model was held at -88 mV and clamped at various voltages from -68 mV through -54 mV in increments of 2 mV following Kang and Kitai (1993b). B: IA, model was held at -90 mV and clamped to various voltages from -60 to 75 mV in 15-mV increments after experiments by Silva et al. (1990). C: Ih, model was held at -57 mV and clamped to various potentials from -141 to -69 mV in 12-mV increments after Mercuri et al. (1995). D: IK, model was held at held at -40 mV and clamped to 20 mV in 10-mV increments after experiments by Silva et al. (1990). E: ICa,N, model was held at -84 mV and clamped to -58, -48, and -38 mV following Kang and Kitai (1993b). Steady-state activation and inactivation curves for IA (<A><AC>p</AC><AC>&cjs1171;</AC></A> and <A><AC>q</AC><AC>&cjs1171;</AC></A>, respectively) and steady-state activation curves for Ih (<A><AC>q</AC><AC>&cjs1171;</AC></A>) and IK (<A><AC>n</AC><AC>&cjs1171;</AC></A>) are displayed in F. Current parameters and equations are in Tables 1-3.

Fits to voltage-clamp data are not the only criteria for formulating ionic current descriptions, which may have to be modified to fit whole cell transmembrane potential data. These additional adjustments are justified considering that the voltage-clamp and the free-running SOP and square-wave recordings were obtained from different cells and in different laboratories and that the voltage-clamp experiments were performed on a particular cell, and there is a considerable variation in the waveshape of the ionic current response from cell to cell.

Calcium currents---ICa

The mathematical equations used in the model for the different calcium currents are given in Table 1, and the voltage dependence of the steady-state activation and inactivation characteristics of these currents are shown graphically in Fig. 3. The reversal potential for the calcium currents has been set to a constant 50 mV by extrapolating current-voltage curves from Kang and Kitai (1993b).


                              
View this table:
[in this window]
[in a new window]
 
Table 1. Calcium currents---ICa



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 3. Steady-state characteristics of calcium currents. d and f represent the activation and inactivation gating variables, respectively, for calcium currents. Steady-state activation (<A><AC>d</AC><AC>&cjs1171;</AC></A>T) and inactivation (<A><AC>f</AC><AC>&cjs1171;</AC></A>T) characteristics of ICa,T were obtained from voltage-clamp studies of Kang and Kitai (1993b). Steady-state activation of ICa,N (<A><AC>d</AC><AC>&cjs1171;</AC></A>N) was obtained by fitting voltage-clamp data from the same source. Steady-state activation of ICa,L (<A><AC>d</AC><AC>&cjs1171;</AC></A>L) is more hyperpolarized than that of ICa,N. Steady-state activation (<A><AC>d</AC><AC>&cjs1171;</AC></A>HVA) and inactivation (<A><AC>f</AC><AC>&cjs1171;</AC></A>HVA) of ICa,HVA, are adapted from Cardozo and Bean (1995) and Kang and Kitai (1993b), respectively.

T-TYPE LVA CALCIUM CURRENT---ICA,T. The T-type calcium current is based on data from Kang and Kitai (1993b). The membrane potential was held at -88 mV and clamped at various voltages from -68 mV through -54 mV in increments of 2 mV. Their experiments show a two-stage inactivation: a quick inactivation to a small current followed by a slow inactivation that maintains the small current for the duration of the voltage clamp (300 ms). Characterization for the long secondary component is based on data from neostriatal neurons (Hoehn et al. 1993), whereas the activation and fast inactivation component was characterized by fitting data from Kang and Kitai (1993b). The voltage-clamp fits to data are shown in Fig. 2A. The current description fits the quick activation and biphasic inactivation well at the expense of matching the peak current values at some clamp potentials.

N-TYPE HVA CALCIUM CURRENT---ICA,N. As modeled, ICa,N begins to activate at -60 mV and continues to activate at -45 mV. The slow inactivation that appears in voltage-clamp records (Kang and Kitai 1993b) is probably calcium-mediated, because no calcium-independent inactivation was detected in the presence of the calcium buffer EGTA and the addition of the calcium buffer EGTA increased the inactivation time constant. We were unable to find specific data regarding the Ca2+-mediated inactivation of ICa,N in DA neurons and hence employed a simple Michaelis-Menten relationship to modify the channel conductance. The parameters of the gating variable equations associated with the equations for ICa,N were determined generally by fitting voltage-clamp data (Fig. 2E) from Kang and Kitai (1993b) in a cell with N-type channels which did not exhibit calcium-dependent inactivation. The cell was held at -84 mV and clamped to voltages of -58, -48, and -38 mV, respectively. Because the characterization of the calcium-dependence of this current was not based on voltage-clamp data, the parameters associated with the steady-state calcium-dependent inactivation function (fCa,N) were loosely set in the fit to voltage-clamp data and determined finally in fits of the complete model to SOP and square-wave data. The specific equations for ICa,N are listed in Table 1.

L-TYPE HVA CALCIUM CURRENT---ICA,L. On the basis of experiments in other cells (Fox et al. 1987a,b; Johnston and Wu 1995), we have assumed that the nifedipine-sensitive, long-lasting calcium current, ICa,L, in DA neurons has activation characteristics that are similar to those of ICa,N. The steady-state half-activation voltage (V0.5) associated with ICa,L was set to be more hyperpolarized compared with the V0.5 value of the ICa,N activation characteristic. This is consistent with observations about the relationship between dihydropyridine- and omega -conotoxin-sensitive currents (Kasai and Neher 1992; Regan 1991) (see Fig. 3).

Recent experiments have shown that calcium chelation with bis-(o-aminophenoxy)-N,N,N',N'-tetraacetic acid (BAPTA) prolongs the plateau of the square-wave oscillation, which is blocked by nifedipine (Ping and Shepard 1997). This is consistent with the hypothesis that the L-type current, like ICa,N, is inactivated by [Ca2+]i. The inactivation of ICa,L has been modeled as calcium dependent (Johnston and Wu 1995). A Michaelis-Menten relationship, (fCa,L in Fig. 4) has been used to characterize the calcium-dependent inactivation of this current, and its parameters have been adjusted to give reasonable fits to SOP and square-wave oscillation data. The equations for ICa,L can be found in Table 1.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 4. Steady-state calcium-dependent characteristics associated with various ionic currents. A: summary of calcium-dependent ionic currents. fCa,L, fCa,N, and pK,Ca are indicated. B: activation of IK,Ca,SK by calcium during slow oscillation potential (SOP). IK,Ca activation variable (pK,Ca) is plotted against the internal calcium concentration, [Ca2+]i. Extent of [Ca2+]i excursion during SOP is indicated by the shaded region. C: inactivation of ICa,L and ICa,N by the calcium-dependent fCa,L and fCa,N, respectively, during square-wave oscillations. Extent of [Ca2+]i excursion during the square-wave oscillations is indicated by the shaded region.

COMBINED HVA CALCIUM CURRENT---ICA,HVA. The remainder of the HVA current has been lumped together as ICa,HVA. This includes the omega -agatoxin-sensitive P-type current and the R-type current, which is not selectively blocked by any known pharmacological agents. The parameters used for voltage-activation (V0.5,act = -10.0 mV, k = 10 mV) are similar to those obtained by Cardozo and Bean (1995), which were -5.0 mV and 13.8, respectively. The parameters for voltage-inactivation are similar to those obtained by Kang and Kitai (1993b), (V0.5,inact = -48 mV, k = 5.0 mV), which were -48 mV and 6.0, respectively.

CALCIUM-ACTIVATED POTASSIUM CURRENT---IK,CA. Dopaminergic neurons are known to contain at least two types of calcium-activated potassium currents, the BK and SK currents (Shepard and Bunney 1988; Silva et al. 1990). The BK, or maxi type channel, has a high channel conductance and voltage- and calcium-dependent activation and is known to be blocked selectively by iberiotoxin and nonselectively blocked by TEA. Specifically, Silva et al. have shown that IK,Ca,BK is blocked by TEA in DA neurons (Silva et al. 1990). The smaller conductance SK channel has a purely calcium-dependent activation and obtains its voltage dependence indirectly via the voltage dependence of calcium entry mechanisms (e.g., calcium channels) (Hille 1992). It is blocked selectively by apamin and is insensitive to TEA. For the purposes of simulation, we have assumed that apamin partially blocks IK,Ca,SK, leaving it at ~20% of its original strength. Further, we have assumed that IK,Ca,SK is activated by intracellular Ca2+ through the Boltzman-type function pK,Ca with a half-activation calcium concentration of KM,K,Ca given in Table 2 and seen in Fig. 4. The calcium half-activation concentration value KM,K,Ca has been set to 190 nM. This is reasonable in view of experimental observations that the maximum calcium excursion of DA neurons during SOP is ~200 nM (Callaway and Wilson 1997). Because IK,Ca,BK is inactivated by TEA application as mentioned earlier, it cannot be essential for the slow underlying oscillations (SOP). In the present modeling study, we have excluded this active spiking mode of oscillation, and therefore we have not included IK,Ca,BK in the model.


                              
View this table:
[in this window]
[in a new window]
 
Table 2. Potassium currents

TRANSIENT OUTWARD CURRENT---IA. The 4-aminopyridine (4-AP)-sensitive current, IA, has been observed in DA neurons and plays a role in the regulation of action-potential frequency by slowing the recovery of membrane potential to baseline levels (Silva et al. 1990). The steady-state activation (p) and inactivation (q) characteristics of IA (Fig. 2F) were determined by fitting published voltage-clamp data obtained from DA neurons of the substantia nigra pars compacta (SNc) of the rat (Silva et al. 1990). The cell was held at -90 mV and clamped to various voltages in the range -60 to 75 mV in 15-mV increments. Equations describing this current are given in Table 2, and Fig. 2B shows the model-generated fits to voltage-clamp data.

HYPERPOLARIZATION-ACTIVATED CATION CURRENT---IH. The hyperpolarization-activated cation current (Ih) has been shown to be important for pace-making activity in thalamic and cortical neurons (McCormick and Pape 1990). In DA neurons, however, Ih has been shown to have a negligible effect on spontaneous firing, resting membrane potential, and normal resting conductance (Mercuri et al. 1995). Our simulations confirm that Ih is not essential for the generation of the SOPs or for setting the normal resting input impedance of the model (see Table 3). However, Ih is the major component of the input impedance of the membrane at more hyperpolarized potentials (about -100 mV).


                              
View this table:
[in this window]
[in a new window]
 
Table 3. Hyperpolarization-activated cation current---Ih

Ih activates with hyperpolarization beyond approximately -60 mV and exhibits no inactivation. Its reversal potential Eh is -35 ± 4 mV (Mercuri et al. 1995), which lies between the reversal potentials for potassium and sodium currents (about -70 and 60 mV, respectively). Ih therefore is modeled as a mixed cation current predominately permeable to sodium and potassium using a previously developed characterization (Demir et al. 1994). Model-generated fits to voltage-clamp data from rat by Mercuri et al. (1995) are shown in Fig. 2C. The cell was held at -57 mV and clamped to various potentials from -141 to -69 mV in 12-mV increments.

DELAYED-RECTIFIER CURRENT---IK. The delayed-rectifier current, IK, was characterized by fitting published voltage-clamp data from rat DA neurons (Silva et al. 1990) in which the cell was held at -40 mV and clamped to 20 mV in 10-mV increments (Fig. 2D). The equations characterizing IK are given in Table 2. As mentioned in the introduction, SOPs are observed after TTX application in the presence IK and enhanced after its blockade by TEA. The significance of this observation is that IK is not necessary for SOP and, unless otherwise noted, is not present in either SOP or square-wave oscillations.

BACKGROUND CURRENT---IB. Although important, there is little quantitative experimental data upon which to base a mathematical description of IB. The linear background current is modeled as a sum of three separate components: a sodium current (IB,Na) a calcium current (IB,Ca), and a potassium component (IB,K) The equations for these linear leak conductances are given in Table 4. The input impedance of the model was obtained by clamping the membrane to -60 mV from a holding potential of -50 mV and observing the evoked transient current. A Delta V of 10 mV was divided by the resultant Delta I to obtain an input impedance of 365 MOmega , which is within the experimentally reported range for these neurons (Grace and Onn 1989; Johnson and North 1992; Yung et al. 1991).


                              
View this table:
[in this window]
[in a new window]
 
Table 4. Background current---IB

PUMP AND EXCHANGER CURRENTS---ICA,P, INA,CA, AND INA,K. The equations (see Table 5) for the sodium-calcium exchanger and the calcium pump are adapted from Canavier et al. (1991) and the weakly voltage-dependent sodium-potassium pump from Lindblad (1996). The value for the half-activation calcium concentration for ICa,pump (KM,CaP) was set to 0.50 µM, which is within the range of experimental observations for other preparations (Lichtman et al. 1981; Michaelis et al. 1983).


                              
View this table:
[in this window]
[in a new window]
 
Table 5. Pumps and exchangers

Fluid compartment model

Calcium dynamics are an important part of DA neuron activity because Ca2+-mediated processes are hypothesized to be involved in the generation of the subthreshold oscillations, and, as a result, the mechanisms of rhythmic firing. The material balance on calcium depends on three processes: entry, extrusion, and buffering of calcium. The differential equation representing the rate of change of cytosolic Ca2+ is
[<A><AC>Ca</AC><AC>·</AC></A><SUP>2+</SUP>]<SUB>i</SUB>=<FR><NU>(<IT>I</IT><SUB><IT>Na,Ca</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>Ca</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>B,Ca</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>Ca,pump</IT></SUB>)</NU><DE><IT>2</IT><IT>Vol</IT><SUB><IT>i</IT></SUB><IT>F</IT></DE></FR><IT>−</IT><IT>n</IT>[<IT>B</IT>]<SUB><IT>i</IT></SUB><IT><A><AC>O</AC><AC>˙</AC></A></IT><SUB><IT>C</IT></SUB> (3)
where Voli is the cytosolic volume in nanoliters, F is Faraday's constant in C/mol, and [B]i is the concentration of the internal calcium buffer in millimolar with n binding sites for Ca2+. nBi gives the total concentration of binding sites for Ca2+. The volume is calculated based on a spherical soma of radius 15 µm for a total intracellular volume of 0.0141 nl.

The buffering of intracellular calcium has been modeled as binding to an intracellular protein such as calmodulin. The formula for the buffer (Robertson et al. 1981) is given by
<IT><A><AC>O</AC><AC>˙</AC></A></IT><SUB><IT>C</IT></SUB><IT>=</IT><IT>k</IT><SUB><IT>U</IT></SUB>[<IT>Ca<SUP>2+</SUP></IT>]<SUB><IT>i</IT></SUB>(<IT>1−</IT><IT>O</IT><SUB><IT>C</IT></SUB>)<IT>−</IT><IT>k</IT><SUB><IT>R</IT></SUB><IT>O</IT><SUB><IT>C</IT></SUB> (4)
where OC is the buffer occupancy or the fraction of sites that already are occupied by Ca2+ and therefore unavailable for binding. All binding sites are considered as a single population and assumed to be independent. This assumption ignores the nonlinearities that occur as a result of multiple cooperative binding sites.


    COMPUTATIONAL ASPECTS
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

The complete system consists of 15 state variables: 3 differential equations describing the time rate of change in membrane voltage (V), internal calcium concentration ([Ca2+]i), and fractional occupancy of the calcium buffer (OC) and 12 other differential equations characterizing voltage-dependent gating variables. Equations 1, 3, and 4, along with the tables for individual membrane currents (Tables 1, 2, 3, 4, and 5); parameter values (Table 6); and initial conditions (Table 7), contain all the information necessary to carry out the simulations presented in this paper. The units used in the model are time in milliseconds (ms), voltage in millivolts (mV), concentration in millimoles/liter (mM), current in nanoamperes (nA), conductance in microsiemens (µS), capacitance in nanofarads (nF), volume in nanoliters (nl), and temperature in degrees Kelvin (°K). All simulations were performed on a Sun Ultra I by forward integration of the coupled system of differential equations using an implicit fifth-order Runge-Kutta method with variable step size designed for stiff systems of differential equations (Hairer and Warner 1990).


                              
View this table:
[in this window]
[in a new window]
 
Table 6. Parameter values


                              
View this table:
[in this window]
[in a new window]
 
Table 7. Initial conditions


    RESULTS
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

In addition to providing a quantitative test that the current descriptions derived from voltage-clamp data are sufficient to reproduce experimental data generated under current-clamp conditions, simulations of the dopaminergic neurons in this study enable us to gain insight into the ionic mechanisms responsible for rhythmic firing as well as to utilize this insight to make experimentally testable predictions. Figure 5 summarizes the model results for the two oscillations of interest and compares model-generated and experimentally observed behavior under different simulated experimental conditions. Because TTX must be present to reveal slow sinusoidal oscillations, INa has not been modeled. TEA application is simulated by letting <A><AC>g</AC><AC>&cjs1171;</AC></A>K = 0. which results in slow oscillatory potentials (SOPs, Fig. 5A). Simulation of the application of apamin blocks most of the outward current IK,Ca,SK and results in square-wave oscillations (Fig. 5B). In either case, the model provides SOPs or square-wave oscillations that agree qualitatively with data from Ping and Shepard (1996).



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 5. Model-generated and experimentally observed oscillations of DA neurons under influence of blocking agents. A: model and data comparison of the SOP. IK shut off to simulate TEA application. B: model and data behavior during oscillations produced by partial blockage of IK,Ca,SK to simulate apamin application. Data from rat scanned and digitized from Fig. 2C of Ping and Shepard (1996).

Slow oscillatory potentials

Figure 6 shows the ionic currents underlying the SOP. Figure 6A shows the SOP (same as Fig. 5A), whereas the intracellular calcium concentration is shown in B. The major currents involved in SOPs are ICa,L and IK,Ca, whereas ICa,N is small relative to ICa,L (Fig. 6C). Note that ICa,T and ICa,HVA (E) are not fully activated because the SOPs occur in a potential range that is outside their range of activity. On the other hand, the calcium inactivation characteristic of ICa,N (Fig. 4A) keeps the magnitude of this current below that of ICa,L in the calcium range of oscillation (nominally, 180-200 nM). The other model currents consist of background and pump currents as seen in Fig. 6D as well as the transient outward current, IA (E). These currents are important in biasing the cell in the appropriate range of membrane potential and cytosolic Ca2+ concentration. The hyperpolarization-activated mixed cation current, Ih, can be eliminated with little impact on the shape and frequency of the SOPs (not shown).



View larger version (26K):
[in this window]
[in a new window]
 
Fig. 6. Underlying currents of SOP. A: membrane potential during SOP. B: [Ca2+]i oscillations conform well to data from Callaway and Wilson (1997). C: ICa,L and IK,Ca are the main ionic currents during the SOP with minor contribution from ICa,N. D: background currents IB, sodium-dependent currents INa,K and INaCa, and the calcium pump, ICa,pump. E: remaining calcium currents ICa,T and ICa,HVA and the transient outward current IA.

Mechanisms underlying SOP

The mechanism underlying the SOP (Fig. 6A) can be understood by considering the calcium- and voltage-dependent processes in the model. Broadly, activation of ICa,L depolarizes the membrane and leads to increased [Ca2+]i (Fig. 6C), which activates IK,Ca,SK (Fig. 4B). Activation of IK,Ca,SK, in turn, hyperpolarizes the membrane, decreases ICa,L, and reduces [Ca2+]i.

Clearly, the calcium transient (Fig. 7C) lags behind the calcium current ICa,L, whereas the potassium current (IK,Ca,SK) is in phase. As IK,Ca,SK increases it limits the membrane depolarization due to ICa,L and institutes a first stage of repolarization (slope L1 in Fig. 7A). The peak of ICa,L is long lasting, and it is only after it begins to decline and concomitantly, IK,Ca,SK peaks, that a second, stronger phase of repolarization is brought about (slope L2 in Fig. 7A).



View larger version (27K):
[in this window]
[in a new window]
 
Fig. 7. Important currents underlying sinusoidal oscillations. Membrane potential during SOP (A) responds weakly (with slope of line L1) at first to the activation of IK,Ca (C). ICa,L (B) remains active, increasing [Ca2+]i (D). Finally, ICa,L regeneratively deactivates, hyperpolarizing the membrane (with slope of line L2).

The model-generated calcium transient during SOP conforms well to experimental observations (Callaway and Wilson 1997). Calcium levels vary within the reported range of 160 and 220 nM, and the calcium oscillations lag the voltage oscillations by ~90°. Nedergaard et al. (1993) have shown that neither Ni2+ nor omega -conotoxin can block the SOP in guinea pig SNc DA neurons, whereas nifedipine application alone is effective. Figure 8 shows that blockade of ICa,T, ICa,N, and ICa,L in the model for the purposes of simulating Ni2+, omega -conotoxin, and nifedipine application, respectively, conforms to these experimental results. The increase in amplitude after blockade of ICa,T is due to the removal of a depolarizing drive and reduction in calcium entry, which allows ICa,L to remain active longer, raising the peak potential.



View larger version (23K):
[in this window]
[in a new window]
 
Fig. 8. Simulation of Ni2+, omega -conotoxin, and nifedipine application. A: control. Application of Ni2+ (B) and omega -conotoxin (C) produce small changes in the character of the SOP, whereas nifedipine application (D) eliminates the SOP.

Square-wave oscillations

The component currents underlying the square-wave oscillations are shown in Fig. 9 in the same order as those for the sinusoidal case (Fig. 6). The major depolarizing currents are again ICa,L and ICa,N; however, the outward K+ current IK,Ca,SK has been blocked. Our model predicts that the mechanism for the square-wave oscillations need not depend on the presence of a distinct hyperpolarization current in this mode of oscillation. Without the activation of IK,Ca,SK by calcium, internal calcium concentration continues to rise and achieves a higher peak level than in sinusoidal oscillations (Fig. 9B), and the frequency of oscillation decreases. ICa,L is inactivated more strongly at the elevated levels of [Ca2+]i. As ICa,L declines due to inactivation, the remaining "residual" currents (IRes, Fig. 10) sum to repolarize the membrane.



View larger version (26K):
[in this window]
[in a new window]
 
Fig. 9. Underlying currents of square-wave oscillations. A: membrane potential during square-wave oscillations. B: [Ca2+]i variation during square-wave oscillations is greater than that seen during SOPs. C: ICa,L is the main ionic current. The residual IK,Ca,SK current and ICa,N are minor contributers in this mode of oscillation. D: calcium and sodium background currents, IB,Ca and IB,Na. E: remaining calcium currents ICa,T and ICa,HVA. F: sodium-dependent pump currents, INa,K and INa,Ca. G: transient outward current (IA), potassium background current (IB,K), and the calcium pump (ICa,P).



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 10. Time-expanded view of square-wave oscillations. A and B: plotted as before and show the membrane potential and cytosolic Ca2+ during square-wave oscillations. C: ICa,L and the sum of the remaining currents (IRes = Itot - ICa,L).

Functional tests

Certain functional tests reveal the robustness of the model and provide explanations of experimental observations. It has been shown that application of TTX alone will result in lower-amplitude and higher-frequency sinusoidal oscillations than those resulting from the coapplication of TTX and TEA (Nedergaard and Greenfield 1992; Kang and Kitai 1993a). Figure 11B is a simulation of membrane potential under TTX application alone. The resulting oscillations in this case are due to interactions between ICa,L and the two outward currents IK, and IK,Ca,SK. Recall from Fig. 7, that due to its calcium-dependent activation, IK,Ca,SK reaches its maximum strength ~100 ms after membrane potential peaks. A voltage-dependent current with relatively fast activation, on the other hand, would lead to a hyperpolarizing drive that coincides with the peak voltage. This drive leads to an earlier onset of hyperpolarization, leading to a lower amplitude and a shorter period.



View larger version (12K):
[in this window]
[in a new window]
 
Fig. 11. Effect of TEA on the model. A: control SOP in the presence of TTX and TEA. B: SOP with bath application of TTX alone results in lower amplitude oscillations.

The mechanisms of the square-wave oscillations also were tested via model simulation and functional tests. Ping and Shepard (1997) have shown that chelation of DA neurons during square-wave oscillation by intracellular injection of the calcium buffer BAPTA prolongs the plateau phase of the oscillation. Because BAPTA is a calcium buffer, we assume that the effect of adding more buffer to that already present in the cytosol is similar to increasing the concentration of the internal calcium buffer [B]i. The value of [B]i used nominally in the model (0.050 mM) was increased to 0.080 mM to simulate calcium chelation. Figure 12 shows the model-generated response to this increase in [B]i. At first, the results seem counterintuitive. An increase in buffer concentration should lead to a decrease in the rate of change of [Ca2+]i, resulting in a smaller excursion in calcium concentration. However, on closer inspection of the figure, we see that the increase in [B]i also leads to a decrease in the rate of ICa,L inactivation, which causes a larger calcium influx and a greater excursion in [Ca2+]i.



View larger version (28K):
[in this window]
[in a new window]
 
Fig. 12. Increased buffering and calcium-mediated currents. A: simulation of chelation by increasing the concentration of the internal calcium buffer results in increased plateau duration as compared with the control case. B: increased buffering delays the rise in [Ca2+]i which allows ICa,L (C) to remain active longer. There is little change in the strength of the other calcium-dependent currents, INa,Ca and ICa,P.

In addition to changes in buffer concentration, the calcium availability can be changed by adjusting other calcium-dependent processes in the model. Two additional currents that have the potential for altering calcium availability in the cytosol are the calcium pump, ICa,pump, and the sodium-calcium exchanger, INa,Ca. Increasing the strength of either calcium extrusion mechanism should decrease cytosolic Ca2+ availability and yield results similar to increasing [B]i. In addition to increasing calcium extrusion, calcium availability also can be restricted. Because we have considered ICa,N nonessential for SOPs, it can be reduced without altering the self-excitability of the model. Figure 13 summarizes the results of these adjustments. The maximal conductances for INa,Ca, ICa,pump, and ICa,N (KNa,Ca, <A><AC>I</AC><AC>&cjs1171;</AC></A>Ca,pump, and <A><AC>g</AC><AC>&cjs1171;</AC></A>Ca,N, respectively) were each changed by 5% in these experiments. We found that INa,Ca alone was able to affect duty cycle, i.e., plateau duration could be changed without a significant change in the period of the oscillation. Figure 13A shows the control square-wave oscillations followed by three different schemes for increasing plateau duration by reducing cytosolic Ca2+. Of the three adjustments, only increasing INa,Ca (Fig. 13D) has the desired effect. Decreasing ICa,N (B) and increasing ICa,pump (C) in fact have very little effect.



View larger version (82K):
[in this window]
[in a new window]
 
Fig. 13. Altering plateau duration. Control square-wave oscillations are shown (A) followed by attempts at increasing plateau duration by decreasing ICa,N (B), increasing ICa,pump (C), and increasing INa,Ca (D). Lightly and darkly shaded regions divide the control period at the minimum, and the dashed line marks the control period from that minimum (P). As the lightly shaded regions show, only increasing INa,Ca had the desired effect of increasing plateau duration. Additionally, INa,Ca adjustment allows for control of plateau duration without significant change in oscillation frequency.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

We have developed a semiquantitative model that can mimic the behavior of DA neurons under a wide range of experimental conditions, including: SOPs in the presence of TTX and TEA; the changes in the period and amplitude of sinusoidal oscillations in response to TEA application; square-wave oscillations in the presence of TTX, TEA and apamin; and the increase in the duration of the plateau phase of square-wave oscillations observed in response to calcium chelation. The key feature of the model is that it can be used to explore the putative mechanisms underlying SOP and square-wave oscillations.

Somatic origin of calcium-dependent oscillations

This model is a representation of the soma only and therefore is a first approximation of the full dynamics in a dopamine neuron, which has a dendritic tree whose branches frequently extend >500 µm from the soma. The data on which the current descriptions are drawn are characteristic of the soma and proximal dendrites because the studies on which the descriptions are based (Cardozo and Bean 1995; Kang and Kitai 1993b; Silva et al. 1990) were conducted in acutely dissociated or thin (150-200 µm) coronal slice preparations. In both of these types of preparations, the dendritic arbor can be truncated. In particular, ICa,T may be localized preferentially on distal dendrites because only 11 of 15 cells in thin coronal slice (Kang and Kitai 1993b) and none in the dissociated preparation (Cardozo and Bean 1995) exhibited this current.

On the other hand, the presence of large calcium currents on the distal dendrites has been demonstrated in thick (400 µM) coronal sections (Dunia et al. 1996). Nonetheless, the assumption that the somatic currents and calcium dynamics drive the SOP is justified by several experimental observations. First, during the SOP the calcium transients in the soma are a substantial fraction of those observed in the dendrites (Callaway and Wilson 1997) despite the large disparity in surface area to volume ratio. Hence there is a significant amount of calcium entry into the soma. Second, the PLD can be evoked by brief depolarizations of the soma, and its rate of activation is strongly dependent on the membrane potential in the soma (Grace and Onn 1989). Finally, one study found that sectioning the distal dendrites of DA neurons had no significant effect on the pacemaker firing frequency (Nedergaard and Greenfield 1992). Assuming the firing frequency is strongly influenced by the SOP, the distal dendrites do not contribute significantly to the SOP as recorded at the soma.

Predictive value of model can aid in experimental design

Although the modeled voltage- and calcium-dependence of ICa,L is consistent with available data, the precise form of the voltage dependence and even the presence of calcium inactivation of this current have not been directly demonstrated in rat dopamine neurons. Hence the model description of this current can be considered a prediction (see following text). In addition, model predictions that INa,Ca controls the duty cycle (Fig. 13D) could be tested if a specific INa,Ca blocker could be applied to the bath.

Homogeneity of DA neurons within and across species

Homogeneity of midbrain DA neurons across species frequently is assumed (Cardozo 1993; Lacey et al. 1989; Sanghera and German 1984; Tepper et al. 1987; Yung et al. 1991) because the PLD, depolarized spike threshold, long-duration action potentials, inward rectification, and other characteristics are observed across species. It has been argued (Kotter and Feizelmeier 1998) that scaling alone does not ensure constant activity patterns between species, but additional compensatory mechanisms (LeMasson et al. 1993; Turrigiano et al. 1994) likely are involved. The model in the current study applies specifically to rat mesencephalic neurons because it is based on data from rat (Cardozo and Bean 1995; Kang and Kitai 1993a,b; Ping and Shepard 1996; Silva et al. 1990).

Within a species, these neurons also exhibit a degree of electrophysiological homogeneity (Cardozo and Bean 1995; Johnson and North 1992; Yung et al. 1991) and similar responses to stimuli (Schultz et al. 1995), although they are located in three adjacent regions: the ventral tegmental area (VTA or A10), the substantia nigra pars compacta (A9), and the retrorubral area (A8). Although we have treated the DA neurons as a homogenous population, some have suggested that there are subpopulations of DA neurons (Nedergaard and Greenfield 1992; Shepard and Bunney 1988), and there is evidence of some heterogeneity in their neurochemistry, pharmacology, and electrophysiology (Roth and Elsworth 1995). In addition, the mix of calcium currents in these cells can vary with age (DeFazio and Walsh 1995).

Does the dihydropiridine-sensitive ICa,L mediate the SOP in all cases?

We have assumed that dihydropyridines block the SOP, based on the following data. Mercuri et al. (1994) reported that 30 µM nifedipine blocked spontaneous pacemaker activity in rats, and by inference, the SOP as well. This is consistent with findings in rat that the calcium conductance mediating the SOP is dihydropyridine sensitive (Ping and Shepard 1997) and with those of Kang and Kitai (1993a,b), who found that the calcium channel blocker Cd abolished the SOP in rat. These findings are also consistent with the results of Nedergaard et al. (1993), who reported that 0.5-20 µM nifedipine abolished the SOP in guinea pig. On the other hand, Fujimura and Matsuda (1989) reported that although the calcium channel blockers Cd2+ and Co2+ as well as Ca-free saline blocked the SOP in guinea pig, 100 µM nifedipine did not. Furthermore Yung et al. (1991) reported that 500 µM Ni2+ abolished the SOP in guinea pig, whereas Nedergaard et al. (1993) reported that 500 µM Ni2+ only slightly attenuated it.

Thus, in guinea pig at least, there are conflicting data. However, the evidence argues for a dominant role for ICa,L. For example, Kang and Kitai hypothesized that ICa,N was responsible for the PLD; however, they did not show that omega -conotoxin abolishes the SOP. Nedergaard et al. (1993) did perform this test in guinea pig neurons and found that 1-10 µM omega -conotoxin did not block the SOP; this casts doubt on the importance of ICa,N in generating these oscillations. The data of Yung et al. (1991) suggest a role for ICa,T because Ni2+ selectively blocks this current (but see Ellinor 1993). Our simulations have shown that it is very unlikely that ICa,T, as characterized by the voltage-clamp data of Kang and Kitai, contributes significantly in the voltage range of the PLD. Similarly, although Kang and Kitai have shown that ICa,N exists in the correct voltage range for the PLD, the simulations show that to generate both SOP and square-wave oscillations under the appropriate conditions, a calcium current with different activation characteristics from ICa,N is necessary. We predict that ICa,L, as we have characterized it, is this missing current. The voltage range of activation of this current is more negative than that typically associated with L-type calcium currents; however, similar L-type Ca2+ channels have been found in several locations, including guinea pig motoneurons (Hsiao et al. 1998), turtle motoneurons (Russo and Hounsgaard 1996), and rat supraoptic neurons (Fisher and Bourque 1996). Further experimental work is required to verify the validity of this prediction and determine the activation characteristics of ICa,L in DA neurons.

It is widely accepted that a persistent, TTX-sensitive sodium current also contributes to the SOP (Grace and Onn 1989; Nedergaard et al. 1993; Ping and Shepard 1996). In the dopamine neurons of the retina, which exhibit similar regular pacemaker activity in a slice preparation (Feigenspan et al. 1998), the sodium current, and not a calcium current, is the dominant pacemaking current. Hence under different circumstances, different currents may contribute to the SOP in differing proportions, perhaps due to self-regulatory compensatory mechanisms that dictate that dopamine cells are intrinsic pacemakers.

Comparison with previous models

Li et al. (1996) developed a model that replicated both NMDA- and apamin-induced burst firing but ignored the SOP and regular pacemaker firing. Their model was conceptual in nature with arbitrary parameters not based on the voltage-clamp data of the particular currents in DA cells. The mechanism that they proposed for NMDA-induced bursting is that same as that proposed by a later, more physiologically based model by Canavier (1999), namely the interaction of the NMDA-induced current (INMDA) and the sodium pump. Regenerative voltage activation of INMDA is postulated to be responsible for the depolarization during a burst, and sodium accumulation due to entry via INMDA is hypothesized to activate the electrogenic sodium pump, which is a net outward current, until a regenerative hyperpolarization due to the voltage-dependent closing of INMDA channels. Sodium is removed during the hyperpolarizing phase, allowing the cycle to begin again. On the other hand, the mechanism proposed by Li et al. (1996) for apamin-induced, calcium-dependent bursting assigns a crucial role to ICa,T rather than ICa,L, which is unlikely in view of the more physiologically based simulations in the current study. In addition to the models described above, Kotter and Feizelmeier (1998) also have modeled DA neurons, but it is difficult to compare the models because that study did not discuss underlying ionic mechanisms, but merely showed that different activity patterns, such as single-spike firing and burst-firing, could arise from scaling parameters to account for different neuronal sizes.

Whereas the model of Canavier (1999) focused solely on sodium dynamics and the current study focuses solely on calcium dynamics, in reality these systems operate in tandem and interact with each other. Further work is required to elucidate this interaction, as apamin is known to facilitate NMDA-induced burst firing (Seutin et al. 1993). As far as the relative contributions of these types of burst-firing to the situation in vivo, it is well known that NMDA promotes burst firing in vivo (Overton and Clark 1992) and that DA neurons receive excitatory amino acid input (Carter 1982; Robledo and Feger 1990; Scarnati et al. 1986). To date, an endogenous substance that modulates IK,Ca,SK in vivo has not been identified. Some have argued that apamin-induced burst firing resembles in vivo firing patterns more closely than burst firing induced by NMDA (Overton and Clark 1997), so a role for IK,Ca,SK may yet be discovered in addition to the well-established role for INMDA.


    APPENDIX A: NULLCLINE ANALYSIS
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Nullcline analysis was an invaluable tool in the development of the model. Repeated, time-consuming simulations were avoided, and considerable insight into the model was gained through nullcline analysis, which is a subset of the more general phase space analysis of nonlinear systems (see Rinzel and Ermentrout 1998).

Analysis of multidimensional systems---application to DA model

A multidimensional system can be analyzed with two-dimensional nullcline methods if it can be reduced to a two-dimensional system. This is done by eliminating the kinetics of all but two state variables and assuming that they can satisfactorily characterize the system. As we have seen, the slow cyclic variation in [Ca2+]i along with membrane potential are essential contributors to oscillation in DA neurons. Therefore we reduce the system to a two-dimensional one and set all state variables to their steady-state values. We obtain the nullclines for the reduced system by setting all state variables to their steady state values and numerically solving the following equations:
<FR><NU>d<IT>V</IT></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=0</IT>

<FR><NU>d[Ca<SUP>2+</SUP>]<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR><IT>=0</IT> (A1)
Additionally, we assumed that a rapid equilibrium is reached between the free calcium ([Ca2+]i) and the buffered calcium (n[B]OC). We then used the ratio of the free calcium to the total calcium to scale the derivative of [Ca2+]i (Chay and Keizer 1983). If we assume that the [Ca2+]i kinetics are slow compared with those of V, we can use the slope of the potential nullcline at the equilibrium point as an indication of the stability of the system. If the slope is negative, we have a stable equilibrium as any perturbations return the trajectory to the equilibrium point. A positive slope indicates an unstable equilibrium point as it repels all points in its vicinity to a limit cycle. The DA neuron nullclines for the SOP and square-wave oscillations along with their trajectories are shown in Fig. A1A. The waveforms generated by the reduced (two dimensional) model (Fig. A1B, - - -) in are compared with those of the complete model (---). In both SOP and square-wave oscillations, the reduced model has a smaller calcium excursion, and as a result, a higher frequency.



View larger version (31K):
[in this window]
[in a new window]
 
Fig. A1. DA neuron nullclines, associated trajectories, and kinetics of the reduced model. Oscillation trajectory in the V-[Ca2+]i plane for SOPs and square-wave oscillations (A). Time course of the potential and calcium oscillations is shown in B. Reduced model (- - -) has a smaller calcium excursion and a higher frequency compared with the complete model (---).

Rough tunings of the model were performed through analysis of the nullclines of the reduced system. First, the slope of the potential nullcline at the equilibrium point of the reduced system was a good predictor of the stability of the complete system. Using this rule, we were able to make large changes in the system equations while maintaining the model in the correct range for oscillation. Also, because [Ca2+]i varies slowly relative to the potential, the trajectory in phase plane remains more or less on the potential nullcline, jumping between stable branches when an unstable region is encountered. This was very helpful in estimating from the nullclines, the extent of [Ca2+]i and V excursion without running the full simulations.


    APPENDIX B: MODEL EQUATIONS AND PARAMETERS
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

The equations describing the dopaminergic neuron model are contained within this appendix. Expressions for the transmembrane currents are given in Tables 1-4. The pump and exchanger currents are given in Table 5. The gating variables in the model are n, dT, fTf, fTs, dN, dL, dHVA, fHVA, p, q1, and q2. The equations for <A><AC>z</AC><AC>&cjs1171;</AC></A>(V) and tau z(V) can be found in Tables 1-3. A listing of model parameter values is given in Table 6, and the initial conditions needed to run the model in both modes of oscillation are given in Table 7 with SOPs and square-wave oscillations in the left and right columns, respectively.


    ACKNOWLEDGMENTS

This work was funded by a Biomedical Engineering Research grant from the Whitaker Foundation and National Institute of Neurological Disorders and Stroke Grant NS-37963.


    FOOTNOTES

Address for reprint requests: C. Canavier, Dept. of Psychology, University of New Orleans, New Orleans, LA 70148.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Received 9 February 1999; accepted in final form 29 June 1999.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
MODEL DEVELOPMENT
MEMBRANE CURRENTS
COMPUTATIONAL ASPECTS
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

0022-3077/99 $5.00 Copyright © 1999 The American Physiological Society