Mathematical model of subthalamic nucleus neuron: characteristic activity patterns and bifurcation analysis

The subthalamic nucleus (STN) has an important role in the pathophysiology of the basal ganglia in Parkinson's disease. The ability of STN cells to generate bursting rhythms under either transient or sustained hyperpolarization may underlie the excessively synchronous beta rhythms observed in Parkinson's disease. In this study, we developed a conductance-based single compartment model of an STN neuron, which is able to generate characteristic activity patterns observed in experiments including hyperpolarization-induced bursts and post-inhibitory rebound bursts. This study focused on the role of three currents in rhythm generation: T-type calcium (CaT) current, L-type calcium (CaL) current, and hyperpolarization-activated cyclic nucleotide-gated (HCN) current. To investigate the effects of these currents in rhythm generation, we performed a bifurcation analysis using slow variables in these currents. Bifurcation analysis showed that the HCN current promotes single-spike activity patterns rather than bursting in agreement with experimental results. It also showed that the CaT current is necessary for characteristic bursting activity patterns. In particular, the CaT current enables STN neurons to generate these activity patterns under hyperpolarizing stimuli. The CaL current enriches and reinforces these characteristic activity patterns. In hyperpolarization-induced bursts or post-inhibitory rebound bursts, the CaL current allows STN neurons to generate long bursting patterns. Thus, bifurcation analysis explained the synergistic interaction of the CaT and CaL currents, which enables STN neurons to respond to hyperpolarizing stimuli in a salient way. The results of this study implicate the importance of CaT and CaL currents in the pathophysiology of the basal ganglia in Parkinson's disease.

The subthalamic nucleus is an important part of the brain due to its involvement in motor behavior and other brain functions.Pathological activity in subthalamic neurons contributes to the motor symptoms of Parkinson's disease.Therefore, it is important to understand the mechanistic basis of electrical activity in subthalamic neurons.These neurons are known to exhibit several distinct patterns of activity: tonic spikes (action potentials), hyperpolarization-induced bursts of spikes, and post-inhibitory rebound bursts of spikes.This study explores how different membrane channels interact with each other to give rise to various activity patterns in a single cell within the framework of a one-compartment (i.e.no spatial extent) model.Bifurcation analysis using slow variables in membrane currents showed how calcium currents can promote different types of bursting, especially in response to an inhibitory input.The study also emphasizes the importance of slow calcium currents in subthalamic cells in Parkinson's disease known for its bursting electrical activity in the subthalamus.
In vitro, the STN neuron fires single spikes in a slow rhythmic manner at 5-20 Hz during the absence of external input, which may underlie the tonic firing patterns observed in vivo in resting animals [Beurrier et al., 1999;Bevan and Wilson, 1999].STN firing frequency increases almost linearly with the magnitude of injected depolarizing current pulses [Hallworth et al., 2003;Wilson et al., 2004].Bevan et al. showed that a majority of STN neurons elicit a calcium-dependent post-inhibitory rebound (PIR) burst of spikes when the neurons are released from inhibitory synaptic input [Bevan et al., 2002a].PIR bursts can be either long or short depending on the level and duration of inhibition received by the neuron.On the other hand, some STN neurons under hyperpolarized conditions switched from a spontaneously discharging single-spike mode to a pure burst-firing mode (consisting of long-lasting bursts of constant duration) or a mixed burst-firing mode (alternating short and long bursts) [Beurrier et al., 1999].It was argued that slow rhythmic bursting results from T-type calcium currents and L-type calcium currents [Beurrier et al., 1999;Bevan and Wilson, 1999].
Since the STN receives inhibitory input from the external globus pallidus (GPe), these characteristic activity patterns may be essential ingredients in the normal and abnormal functioning of the system.For example, STN neurons are able to transform sustained inhibitory synaptic input into rhythmic bursts of spikes.Hence, their resting activities can transition from tonic spiking discharge to burst-firing patterns.Therefore, hyperpolarization-induced bursting in STN neurons may play crucial roles in the generation of excessively synchronized rhythmic bursting patterns observed in parkinsonian BG (see references above).
The goal of this study is to develop a relatively simple STN neuron model which exhibits several activity patterns as described above, and to understand the dynamic mechanisms of these patterns by studying the interactions of membrane currents and the bifurcation diagrams underlying transitions between different activity modes.A model, which is able to generate all types of activity, was developed earlier [Gillies and Willshaw, 2006], but it is a complicated multi-compartmental model where the interaction of compartments appears to be essential for its dynamical regimes.While this model is able to generate the activity patterns described above, mathematical analysis of the underlying mechanisms is very challenging due to the complexity of the multi-compartmental model.In addition, the model's complexity makes it hard to use in the construction of a large model network, such as the cortex-BG network to study the mathematical mechanisms of synchronized beta rhythms in networks.
In our current study, we developed a conductance-based single compartment model of the STN neuron and showed how this model captures the characteristic activity patterns, especially hyperpolarization-induced bursting and post-inhibitory rebound (PIR) bursting.Using this model, we performed a bifurcation analysis to study the roles and effects of three currents (T-type calcium current, L-type calcium current, and hyperpolarizationactivated cyclic nucleotide-gated current) in the rhythm generation mechanisms under inhibition [Beurrier et al., 1999;Atherton et al., 2010].

Mathematical model
A conductance-based single compartment model of the STN neuron includes spikegenerating potassium and sodium currents (IK and INa), a leak current (IL), a persistent sodium current (INaP), a calcium dependent potassium current (IAHP), hyperpolarizationactivated cyclic nucleotide-gated (HCN) current (IHCN), an A-type potassium current (IA), a T-type low-threshold calcium (CaT) current (ICaT), and an L-type high-threshold calcium (CaL) current (ICaL).For basic currents (spike-generating currents, leak currents, and AHP currents), we used the equations in the Terman model [Terman et al., 2002;Best et al., 2007].The forms of A-type calcium current, CaT current and CaL current were adopted from the Hahn and McIntyre model [Hahn and McIntyre, 2010].The form of the HCN current was based on the Gillies and Willshaw model [Gillies and Willshaw, 2006].The dynamics of the membrane potential (V) is described by the following differential equations: ,  ∈ {,  7 } (3) where ).The units for ionic currents are mA/cm 2 .In the first equation,  '--. is the baseline external input, and  '--represents an injected applied current.In the last equation, [Ca] is the calcium concentration in mM, F is the Faraday's constant and KCa = 0.2/ms is the calcium pump rate.Voltage dependent activation and inactivation steady states and time constants are given as follows:

24
,  ∈ {, ℎ, , , , , , , ,  4 } (5) ∈ {, ℎ, , , , , , , ,  4 ,  7 } (8) @ () =  .,@+  4,@ GexpH 4,@ +  4,@ K + expH 7,@ +  7,@ KL 24 (9) The values of maximal conductances are as follows:  [Hahn and McIntyre, 2010] and the Gillies and Willshaw model [Gillies and Willshaw, 2006], then tuned to capture the characteristic activity patterns, especially hyperpolarization-induced bursts and PIR bursts.In the following sections, we will focus on CaT, CaL, and HCN currents and study the effects and interactions of these currents on hyperpolarization-induced rhythms.In the next section, we begin with spontaneous tonic firing activity.Although spontaneous tonic spiking is not induced by hyperpolarization, the results in that section will be used to explain the effect and interaction of CaT and CaL currents on hyperpolarization-induced rhythms.Note that in the rest of the manuscript, we will omit the unit of each parameter for the sake of brevity and simplicity.Bifurcation parameters in the following sections are formed by multiplying maximal conductances and dimensionless gating variables, thus they follow the same units as the maximal conductances used therein.
As observed in experiments and discussed in other studies [Hallworth et al., 2003;Wilson et al., 2004], firing frequency increases almost linearly as the magnitude of injected depolarizing current increases (Fig. 1A) or as  +', increases (Fig. 1B).There is only a slight change in frequency when  *+& increases (Fig. 1C).In the case of CaL current variation, however, the frequency shows an abrupt jump around  +'$ = 17 while the frequencies are almost constants for small or large values of  +'$ (Fig. 1D).In Sec 3.2, we will explore how the CaL current and other currents interact to generate these firing patterns through a bifurcation analysis.Note that, when we increased the timescale of the HCN gating variable  dynamics by 50% (that is, the value of  @ was multiplied by 1.5.Note that  .,@= 0.) in Fig. 1D, there is almost no change in frequency for a wide range of  +'$ values.In Sec 3.2, we will discuss this observation in detail.
The mechanism underlying the effect of applied current on spontaneous tonic firing activity has been extensively studied in different types of neuron models [Terman, 1992;Ermentrout and Terman, 2010;Park and Rubin, 2013].The response of these neuron models and the model considered here to the constant applied current are similar and have been considered in the studies mentioned above.Therefore, we will focus on the effect of the three currents considered here (CaT, CaL, and HCN currents) on spontaneous tonic firing activity.In the following sections, we study how the role and mechanisms of these currents interplay in each characteristic firing pattern using bifurcation analysis.We will utilize the separation of timescales technique and slow variables to construct bifurcation diagrams (the fast-slow analysis).

Effect of CaT and HCN currents on spontaneous tonic firing rhythm. Bifurcation analysis.
To study the effects of the CaT current and the HCN current on spontaneous tonic firing rhythms, we used a fast-slow analysis with  +'$ fixed at 5. In the standard fast-slow analysis, we consider slow variables as bifurcation parameters and derive a bifurcation diagram of the fast subsystem.During spontaneous tonic firing activity in the model, the level of [] is very low and changes slowly, which is due to the slow timescale of [] and low spiking frequency.Hence, we may disregard the effect of [] and the CaL current in this case.Note that the gating variables of the CaT current and the HCN current do not depend on [].On the other hand, it is known that the HCN current promotes a single-spike activity [Atherton et al., 2010].For this reason, we used the following slow variables ℎ =  *+& *  and  =  +', *  for the analysis, where f and q are the gating variables in HCN and CaT currents, respectively.In summary, f and q (or ℎ and ) are slow variables and the remaining variables are then fast variables.The system of governing equations of the fast variables forms a fast subsystem.
Figure 2A shows one example of a bifurcation diagram of the fast subsystem, which was projected onto the (, )-space.Here we treated  as a bifurcation parameter and fixed ℎ at 0.01.The set of fixed points of the fast subsystem forms an S-shaped curve, , in the (, )-space and this structure persists over a certain range of ℎ values of interest.The lower branch of  at the lower left corner (black solid) consists of stable fixed points and the middle branch of  unstable saddle points (black dashed).The lower branch and middle branch coalesce at a fold bifurcation, which we call the right knee (RK) of .Similarly, this middle branch turns around at another fold bifurcation point.We call this upper fold bifurcation point the left knee (LK) of .As  increases from LK, the fast subsystem undergoes an Andronov-Hopf (AH) bifurcation, at a value of  that we denote by  )* , and above this, the fixed points become stable.A family of stable periodic orbits (P) emerges from  at the AH-point.The two black curves from the AH-point show the minimum and maximum V along the family of periodic orbits.Finally, the family of stable periodic orbits (P) terminates in a saddle-node on invariant circle (SNIC) bifurcation.Figure 2A also shows the projection of the spontaneous tonic firing solution (blue) of the full model onto the bifurcation diagram when  +', = 20, along with the corresponding nullcline (green).To the right (above) of the -nullcline,  B < 0; hence,  decreases over that region as a dynamic variable.Similarly,  B > 0 to the left (below) of the nullcline.Note that -nullcline lies above the lower branch of  in the diagram.As will be explained later in this section, two characteristics of the bifurcation diagram --1) the family of stable periodic orbits terminates at the SNIC and 2) the -nullcline lies above the lower branch of --are crucial for the generation of tonic spiking solutions.For this reason, this type of bifurcation has been frequently associated with tonic spiking solutions.As will be explained in Section 4 (cf.Figures 4 and 5), if the stable periodic orbits terminate in a homoclinic orbit, then a square-wave bursting in the neural model is more likely to occur.[Rinzel, 1987;Best et al., 2005;Butera et al., 2005;Bertram and Rubin, 2017] Now treating  and ℎ as bifurcation parameters, a bifurcation diagram in Fig. 2A forms a bifurcation surface in the (, ℎ, )-space (Fig. 2B).The black surfaces on the top and bottom of the right side of the figure are surfaces formed by maximum and minimum values of V along the families of stable periodic orbits (cf.Fig. 2A).Let S denote the surface formed by S-shaped curve  from the Fig. 2A.Part of the surface S is also shown between the surfaces of stable periodic orbits.More specifically, the black surface that crosses the figure horizontally in the middle is the surface of unstable fixed points.The folded surface corresponds to the lower and middle branches of S. This surface is folded at the line of RK points, separating stable and unstable points.Figure 2B also shows the projection of the spontaneous tonic firing solution when  +', = 20 (blue), the same solution as used in Fig. 2A.The  -nullsurface is also shown in green.This nullsurface divides the phase space into two parts.In the region below the -nullsurface,  B > 0, hence  increases.On the other hand, in the region above the surface,  B < 0 and  decreases.
The projection of tonic spiking solution sweeps both regions.Figure 2B shows that the projected trajectory jumps down to the lower part of S. Consider a projected trajectory that lies in a small neighborhood of the lower part of S (surface of stable fixed points).Because the lower part of S lies below the -nullsurface ( B > 0),  increases along the surface.In this case, ℎ also increases and as a result, the trajectory traverses the lower part of S. Once the trajectory reaches the line of RK points (the boundary between the lower and middle parts of S), then the trajectory jumps up into the regime of stable periodic orbits of fast subsystem.While staying in the regime of stable periodic orbits, both  and ℎ decrease, hence the trajectory moves away from the regime of stable periodic orbits, and then it jumps down to the lower part of S, which completes one cycle of action potential.
We note that the existence of spontaneous tonic firing solutions depends on the following two facts: 1) the family of stable periodic orbits of the fast subsystem terminates in a saddle-node on invariant circle (SNIC) and 2) -nullsurface lies above the lower part of S. If  +', is decreased, then the -nullsurface moves downward and intersects with the lower part of S while the family of stable periodic orbits still terminate in a SNIC.The resulting intersection curve is the curve of globally stable fixed points.Hence, the trajectory approaches this curve and cannot enter the regime of stable periodic orbits.As a result, there are no spontaneous tonic firing solutions for small  +', values.
Two-parameter bifurcation diagram is given in Fig. 2C, which shows RK line (black slant line) with the projection of four spiking solutions for  +', = 15 (black), 20 (blue), 25 (red), and 40 (magenta).At upper right turning point, trajectory turns counterclockwise.Figure 2C shows that the projection of tonic spiking solution becomes flat and moves rightward as  +', increases.Based on these two observations, we can provide a heuristic explanation for why the frequency of spiking solution increases as  +', increases as follows.First, note that  increases faster for larger  +', values because  +', controls the dynamics of  =  +', * .On the other hand, we may assume that the dynamics of ℎ is similar in these four examples because  *+& is fixed for all cases ( *+& = 2).Thus, the period of spontaneous tonic firing solution may be estimated by the range of ℎ in the trajectory.Thus, we can expect that the frequency tends to increase and to level off eventually as  +', increases.Second, recall that, once the trajectory jumps down, then it approaches the lower part of S, moves slowly along it, and then eventually crosses the RK line to jump up.In the two-parameter bifurcation diagram, this corresponds to the portion of the trajectory from lower turning point to the upper turning point.In fact, trajectory spends most of time near the lower part of S. Now, as  +', increases, the projected trajectory moves rightward, hence it spends less time near the lower part of S due to the proximity to the RK line.In summary, as  +', increases, trajectory spends less time near the lower part of S with slower dynamics, which results in the increase of the spiking frequency.
Now, we investigate why the projection of tonic spiking solution moves rightward as  +', increases.First, recall that (, ℎ, ) phase space is divided into two subregions by  -nullsurface (Fig. 2B).The place where spiking solution is found seems to be determined by the balance between the times that trajectory spends in these two regions while traversing the phase space.Here we note that the position of -nullsurface in Fig. 2B depends on  +', value.In fact, numerical simulation shows that -nullsurface moves up as  +', value increases.Thus, the increase of  +', value affects the position of spontaneous tonic firing solution near RK.
Since  =  +', * , the position of the spiking solution might be determined by the dynamics of slow variable  in the active phase.Averaging method is frequently used to obtain the reduced autonomous equation that governs the evolution of the slow variable in the active phase.Formally, if = (, ) is the system for the evolution of slow variable  and (, ) is of period T, then the associated autonomous averaged system in the active phase is given by

∶= 𝑓(𝑥̅ ) bbbbbb
, which describes the dynamics of slow variable in the active phase.Now, 3 & (") , hence the position of the spiking solution might be determined by the averaged value of  < ().In this case, since the time constant  C () also depends on the voltage, we used a weighted averaging method for  < () in the regime of stable periodic orbits.More precisely, we computed the following averaged value of  < () in a periodic regime as where T is the period of a stable periodic orbit.That is, for a fixed ℎ value, we draw a bifurcation diagram and find  value for SNIC, say  D&E+ .From  D&E+ , we consider 30 points of  with step size 0.001.Here, each fixed  value corresponds to a periodic orbit.Now we computed averaged  < bbbb over this periodic orbit and Fig. 2D shows the averaged  ( bbbbb =  +', *  < bbbb), when ℎ = 0.01.The results for other ℎ values show qualitatively similar patterns.We checked three  +', values, 20 (blue), 30 (red), and 40 (magenta).Diagonal solid line is the line of identity.Since  bbbbb denotes an equilibrium value of  over a spiking solution, if  bbbbb line lies above (below, resp.) the line of identity, then  is forced to increase (decrease, resp.).Thus, the intersection between the line of identity and  bbbbb curve denotes the  value where spontaneous tonic firing solution tends to reside.As shown in Fig. 2D, the intersection point moves rightward as  +', increases and this result explains why the projected trajectory in Fig. 2C moves rightward as  +', increases.We also note that when  +', = 40 in Fig. 2D, the intersection point is around 0.125 and Fig. 2C shows that the value lies on the right side of RK line.This suggests that the spontaneous tonic firing solution in this case may lie inside the regime of stable periodic orbits and the trajectory spends almost no time on the lower surface of S which results in a higher frequency.

Effect of the CaL current on spontaneous tonic firing. Bifurcation analysis.
In this section, we present the effect of the CaL current on spontaneous tonic firing by varying  +'$ values (default value is 5) with fixed  +', and  *+& values.We tested six different  +'$ values (5, 10, 15, 20, 30, and 40).Over these  +'$ values, we found that the model yields low frequency spiking solutions when  +'$ ≤ 15 and high frequency spiking solutions when  +'$ ≥ 20 (Fig. 1D).As seen in Fig. 1D, the frequency of the spiking solution was not affected by the HCN current even if the timescale of the dynamics of the HCN current is substantially increased.That is, the timescale of kinetics of the HCN current may not play a significant role when the CaL current is not negligible anymore.For this reason, we utilized the following slow variables for bifurcation analysis:  =  +', * ,  =  +'$ *  4 *  7 , and [] , where  4 and  7 are slowly activating/de-activating gating variables in the CaL current, and [] is the calcium concentration.We would like to note that the numerical simulation shows that average calcium level [] is around 0.05 in low frequency spiking solutions and is around 0.18 in high frequency spiking solutions.(red).Green curves are -nullcline when  +', = 20.There are several things to note on these bifurcation diagrams.First, we see that the lower and middle parts of the S-shaped curve of fixed points () remain almost the same over various  values.Especially the dependence of  value of the right knee (RK) of  ( F% ) on  values is negligible.This is clearly shown in Fig. 3C, where the vertical lines denote the RK lines for [] = 0.05 (black) and [] = 0.18 (blue).Second, there are two different termination mechanisms of stable periodic orbits depending on the  values.As shown in the previous section, the family of stable periodic orbits emanates from  at the AH-point.For smaller  values such as 2, stable periodic orbits terminate in a saddle-node on an invariant circle (SNIC) bifurcation.On the other hand, for larger  values such as 6 and 10, the stable periodic orbits turn around at saddle-node bifurcation of periodic orbits (SNPO) to become unstable periodic orbits.The third thing to note is the relative position of -nullcline with respect to , especially with respect to the RK of .As seen in the figure, -nullcline lies above the RK for smaller [] values (Fig. 3A).But -nullcline intersects at the middle and lower parts of  for larger [] values (Fig. 3B).
The former case is similar to the one in the previous section, which implies that there will be a spontaneous tonic spiking solution near the RK for smaller [] values.For larger [] values, on the other hand, a different mechanism comes in.Note that the intersection between -nullcline and the S-shaped curve on the lower branch of stable fixed points is a globally stable fixed point.Hence, if a trajectory jumps down to the lower branch of , then it will approach this globally stable fixed point along the lower branch and remain there.Thus, there is no spiking solution near the RK.The other possible way to obtain a spiking solution is inside the regime of stable periodic orbits if the averaged  value ( bbbbb ) is somewhere between the SNPO and the RK.Since the trajectory does not lie on the lower stable fixed points and does not stay in the periodic orbit, the frequency is higher than ones near RK (Fig. 3A).
Figure 3C shows a two-parameter bifurcation diagram with the projection of full model spiking solutions.The two vertical lines in the figure are the lines of the RK for [] = 0.05 (black) and 0.18 (blue).As shown in Fig. 3A-B, these lines are almost independent of  values.Dotted lines that emanate from the RK lines are SNPO lines.To obtain these SNPO lines, we checked the  values of SNPO for  = 1, 2, 3, …,10, and then connected these points.Thus, the SNPO curves are not smooth enough at some places.Black horizontal lines on the black SNPO line denote the projection of full model spiking solutions for  +'$ = 5, 10, and 15 (from bottom to top).The short horizontal lines at the upper left part denote the projection of full model spiking solutions for  +'$ = 20 (black), 30 (blue), and 40 (red).We can see that spiking solution resides near the RK for small  +'$ values and occurs away from the RK line for large  +'$ values.This suggests that for small  +'$ values, the trajectory jumps down to the lower branches of the stable fixed points and approaches to RK line while for large  +'$ , the trajectory does not jump down to the lower branches of the stable fixed points.
Figure 3D shows a bifurcation surface in (, , )-space when [] = 0.18 with the projection of tonic spiking solutions for  +'$ = 20 (black), 30 (blue), and 40 (red).Minimum and maximum values of V along the family of periodic orbits form black surfaces in the figure, which are folded at the SNPO lines.Inner parts correspond to stable periodic orbits and outer parts to unstable periodic orbits.The S-shaped surface of fixed points (S) is shown in cyan.Figure 3D shows that spiking solutions exist inside the regime of stable periodic orbits.Since slow variables (, ) are not slow enough, spiking solutions are not confined inside the region between two surfaces of stable periodic orbits.As  +'$ decreases, spiking solution approaches the SNPO line.
Figure 3E shows the averaged  value ( bbbbb ) for various [𝐶𝑎] and  values.Recall that the intersection between the averaged  bbbbb and the line of identity is where spiking solution tends to reside.The upper two curves are for  = 2 and 4 when [] = 0.05.These cases are similar to those in previous sections so that a spiking solution tends to reside near the RK.The lower two horizontal curves are for  = 8 and 10 when [] = 0.18.This result shows that spiking solution tends to reside inside a periodic regime.Since [] increases while voltage goes up, the overall level of [] depends on the frequency of spiking solution.When spiking solution is near the RK, the spiking frequency is low, hence, the overall level of [] is low, too.Therefore, when [] = 0.05 (Fig. 3A),  bbbbb values at upper right corner have values close to what is expected (Fig. 3E).On the other hand, when [] =0.18 (Fig. 3B),  bbbbb have values at lower left corner (Fig. 3E).
In summary, when  +'$ is relatively small, the model yields a spontaneous tonic spiking solution, which is similar to those shown in previous section and is facilitated by the CaT current.But, when  +'$ is sufficiently large, then the CaL current pushes the trajectory into the regime of stable periodic orbits, hence there is an abrupt jump in the frequency during this transition and we obtain higher frequency spiking solutions.In terms of a bifurcation diagram, the availability of the CaL current pushes a bifurcation diagram upward, hence the RK of S-shaped curve  lies above the -nullcline.This change does not allow for a spontaneous tonic spiking solution near the RK but creates a solution inside the regime of stable periodic orbits.

Hyperpolarization-induced bursting rhythms
Experimental and computational studies point to the importance of CaT and CaL currents for the generation of hyperpolarization-induced bursting rhythms in STN.Beurrier et al. [Beurrier et al., 1999] showed that some STN neurons can switch from the spontaneous tonic firing to slow bursting rhythms or mixed burst-firing patterns under a sustained hyperpolarizing current application.They argued that CaT and CaL currents underlie the generation of the slow rhythmic bursting.Gillies and Willshaw [Gillies and Willshaw, 2006] showed that their multi-compartment model generates a slow rhythmic bursting in the presence of a uniform reduction in the Ca-dependent SK conductance (simulating the application of apamin) and constant hyperpolarizing current injection.They argued that the interaction of CaT and CaL currents determines the presence and nature of the rhythmic bursting.They also argued that a sufficiently strong CaT current was necessary for the generation of individual bursts.Our model is also able to generate bursting rhythms under a sustained hyperpolarization.Similar to the model in [Gillies and Willshaw, 2006], it was necessary to reduce the AHP current conductance and increase the CaT current conductance to generate bursting rhythms.In this section, we study the effect of CaT and CaL currents as well as the HCN current on hyperpolarization-induced bursting rhythms via bifurcation analysis.

Effect of CaT and HCN currents on hyperpolarization-induced bursting rhythms. Bifurcation analysis.
First, we will explore how CaT and HCN currents in an STN neuron under hyperpolarization generate bursting rhythms without the CaL current.The default parameter values in this case are  +', = 25,  )*( = 0.2,  +'$ = 0,  *+& = 2, and  '--= −16.Note that  '--is a large (in terms of the magnitude) negative number to generate hyperpolarization-induced bursting rhythms.
Figure 4A-B show bursting rhythms for different values of  +', with fixed  *+& (Fig. 4A) and for different values of  *+& with fixed  +', (Fig. 4B).As  +', increases, the period decreases while burst duration and number of spikes within a burst increase (Fig. 4C).Thus, the inter-burst interval decreases too.On the other hand, as  *+& increases, we observed that period, burst duration, and number of spikes within a burst decrease at the same time (Fig. 4C).Please, note that the checkerboard pattern shown in the middle figure (burst duration) is due to the typical spike-adding procedure which is frequently found in bursting regime with a small number of spikes within a burst.
To explore the underlying mechanisms of these results, we performed a bifurcation analysis using  and ℎ (as defined in the previous section:  =  +', *  and ℎ =  *+& * ). Figure 4D shows bifurcation diagrams of the fast subsystem projected onto (, )-space with a bifurcation parameter  for ℎ = 0.1 (black) and 0.2 (red).
Green curve is the -nullcline for  +', = 25.While bifurcation structures are qualitatively similar to those shown in previous section, there are some important characteristic differences compared to those in the previous section.First, the family of stable periodic orbits lies above -nullcline.This holds true for all reasonable ℎ values.Hence, when the projection of full model solution is in a bursting mode, or in other words, when the projected trajectory is inside the regime of stable periodic orbits, the trajectory moves leftward and eventually jumps down to the lower branch of the S-shaped curve of the stable fixed points () at the homoclinic (HC) point.Second, bifurcation diagrams show that middle and lower parts of  move leftward (to the lower values of ) as ℎ increases.When ℎ is small, there is an intersection point between -nullcline and the lower branch of , which is globally stable.As ℎ increases, this intersection point approaches RK of  .Hence, when ℎ is small, if the trajectory is in the small neighborhood of the lower branch of , then it moves along the lower branch of  to approach this intersection point.On the other hand, while approaching, ℎ value increases, which results in the loss of this globally stable fixed point.Then, the trajectory is able to jump up into the stable periodic orbit regime and spiking within a burst begins.
Figure 4E shows bifurcation surfaces in (, ℎ, )-space with bifurcation parameters  and ℎ.Projection of bursting solution for  +', = 25 and  *+& = 2 is also shown (red).The  -nullsurface is omitted and the S-shaped surface of fixed points (S) is shown in cyan for the clarity of the figure.The projection of bursting solution is not confined inside the regime of stable periodic orbits since slow variables (, ℎ) are not sufficiently slow.If we make  and ℎ slower, then we can obtain bursting solutions with a longer burst duration and a larger number of spikes, which are confined inside the regime of stable periodic orbits.However, current bifurcation diagrams are sufficient to analyze hyperpolarization-induced bursting mechanisms.Figure 4E shows that the projected trajectory moves along the lower part of S; both  and ℎ values increase.Once the trajectory crosses the RK line, it jumps up into the regime of stable periodic solutions.
Since the family of stable periodic orbits lies above the  -nullsurface,  and ℎ keep decreasing while the trajectory is inside the stable periodic orbit regime.The trajectory eventually jumps down to the lower part of S, which completes one cycle of a bursting solution.
Figure 4F shows a two-parameter bifurcation diagram with projection of full model bursting solutions for different values of  +', with fixed ℎ (all values are the same as in the Fig. 4A with time-series).Recall that both  and ℎ increase over the silent phase of bursting solution and this corresponds to the almost straight, increasing part of the projected trajectory.Once it passes the RK line, the trajectory turns around counterclockwise to start spiking.Wiggles of the projected trajectory near its lower left part correspond to the active phase of a burst.If the trajectory crosses the HC line, it jumps down to the lower branch of stable fixed points and the active phase of a burst is terminated.Note that the active phases of burst in all three trajectories terminate at similar places (cf.Fig. 4F).As  +', increases, the projected trajectory is stretched horizontally and it moves down rightward, hence the range of  increases, whereas the range of ℎ decreases.We may assume that ℎ evolves on similar timescales in all three bursting solutions when  *+& is fixed.Thus, the period and the interburst interval of the bursting solution can be estimated by the range of ℎ.More specifically, they can be estimated by the maximum value of ℎ where the trajectory turns around to begin the active phase of a burst in the figure.Since this maximum value of ℎ decreases as  +', increases, we expect that the period and the interburst interval also decrease as  +', increases.
A similar argument can be applied to the bifurcation diagram in Fig. 4G, which shows the projection of a full model bursting solutions for different values of  *+& and fixed  +', (all values are the same as in Fig. 4B with time-series).In this case, as  *+& increases, the projected trajectory is stretched vertically and moves up leftward.Thus, the range of  decreases whereas the range of ℎ increases.When  +', is fixed, the interburst interval and the period can be estimated by the range of  since we may assume that  evolves on similar timescales in all three bursting solutions.Therefore, the interburst interval and the period decrease as ℎ increases because the range of  decreases as ℎ increases.
The burst duration and the number of spikes within a burst are determined by the proximity of the trajectory to the AH line.Note that proximity to the AH line means there is a higher chance to fire because all trajectories terminate at similar places for different values of  +', and  *+& (Fig. 4D-G).As a result, the burst duration, and the number of spikes within a burst increase in both cases when trajectory moves rightward or more precisely, as  +', increases with fixed ℎ and as ℎ decreases with fixed  +', .
Figure 4F,G also tells us about the role of the HCN current.Without the HCN current (ℎ = 0), the CaT current driven bursting solution is not possible.Recall that the active phase of a burst is initiated when the trajectory crosses the RK line.Without the HCN current, the projection of solution trajectory approaches the globally stable fixed point on the lower branch of the fixed points (Fig. 4D), hence the trajectory cannot cross the RK line to initiate a burst.In two-parameter bifurcation diagrams, the lack of the HCN current means that the projected trajectory lies on the horizontal axis (Fig. 4F-G).Similarly, we can argue that small HCN current means a very long interburst interval.Also, note that the number of spikes within a burst is determined by the proximity to the AH line.Consequently, due to the shape of the AH line, the number of spikes within a burst will be capped.Large HCN currents would result in high frequency spiking since the interburst interval and the number of spikes within a burst decrease at the same time.This result concurs with the experimental observations showing that HCN channels promote single-spike activity rather than bursting rhythms [Atherton et al., 2010].

Effect of the CaL current on hyperpolarization-induced bursting rhythms. Bifurcation analysis.
In this section, we studied the effect of the CaL current on the hyperpolarization-induced bursting rhythms.Default parameter values are  +', = 25,  )*( = 0.2,  +'$ = 15,  '--= −22,  *+& = 2. Figure 5A shows voltage profiles of bursting solutions for different values of  +'$ .Here, we found that burst duration increases substantially as  +'$ increases, while interburst interval increases slightly.These results are summarized in Fig. 5B, which shows period (solid line), interburst interval (dashed line), and burst duration (dotted line) as a function of  +'$ .While both burst duration and interburst interval increase with  +'$ , the increase of burst duration is more significant and so is the period.As compared to CaT current bursting solutions considered in the previous section, the burst duration and the number of spikes within a burst are substantially increased when the CaL current is turned on.Figure 5C,D shows period and burst duration over two-parameter space for  +'$ = 5, 15, and 25.As shown in Fig. 5C, the region for large periods increases as  +'$ increases.
To explore the underlying mechanisms that result in these differences, we performed bifurcation analysis using  =  +'$ *  4 *  7 ,  =  +', *  , and [].We did not choose ℎ as a slow variable although ℎ is slow over silent phase of bursting solution.In fact, if we increase the timescale of the HCN current, we still obtain qualitatively similar results.This may be because the HCN current is mostly involved in the transition from the silent phase to the active phase of bursting solution while the CaL current affects the active phase of the bursting solution more significantly.Figure 5E shows two exemplary bifurcation diagrams with bifurcation parameter  for [] = 0.4 and  = 5 (dotted), and [] = 0.4 and  = 10 (solid).The green curve is a -nullcline for  +', = 25.Similar to the bifurcation diagrams shown in Fig. 3, we observed that 1) the lower and middle parts of S-shaped curve of fixed points () remain almost the same over various  values, 2) the AH-point moves leftward as  increases, and 3) the stable periodic orbits turn around at saddle-node bifurcation of periodic orbits (SNPO) to become unstable periodic orbits.On the other hand, there are two important differences to note between these diagrams and the ones in Fig. 3. First difference is the relative position of -nullcline with respect to the stable periodic orbits.Bifurcation diagrams in this section show that the branches of stable periodic orbits lie above -nullcline for small  values.In this case, bursting solution keeps moving leftward over the active phase of a bursting solution.For large  values, the stable periodic orbits intersect with -nullcline and the averaged  ( bbbbb ) will have a value close to  value at SNPO (cf.Fig. 3).The second difference is the proximity of SNPO to -nullcline for large  values.This proximity forces the bursting solution to slow down while approaching SNPO over the active phase of bursting solution.These two facts imply that bursting solution slows down while approaching  bbbbb value near SNPO and contribute to longer active phase of bursting for large  values.
Figure 5F shows SNPO surface (left), the AH surface (middle slant), and the RK surface (right) in (, , [])-space.The same figure also shows projections of the three bursting solutions for different values of  +'$ .Note that large  +'$ values mean there is an elevated range of  values.Over the silent phase of a burst,  and  increase while [] decreases.In this model, [] decreases sufficiently fast so that the RK line when [] = 0 roughly determines when the cell crosses the RK line and enters the regime of stable periodic orbits.Similar to the bifurcation diagrams in Section 3.2, we see that the RK surface is almost vertical with fixed [] value, in other words, almost independent of  values.This fact implies that the bursting solution is facilitated by the CaT current.Once the cell enters the regime of stable periodic orbits,  and  begin to decrease while [] begins to increase initially and then decrease over an active phase of burst.Spiking within a burst is terminated once trajectory crosses the surface of SNPO.Since burst duration is determined by the distance between the RK surface and the SNPO surface, burst duration increases as  +'$ increases.Especially because, as evidenced by Fig. 5E, bursting solution spends longer time near the SNPO before jumping down for larger  values, thus there is even longer bursting duration and large number of spikes within a burst for large  values.

Post-inhibitory rebound burst
In this section we study calcium-dependent post-inhibitory rebound (PIR) bursts of spikes when the model neuron is released from application of inhibitory current.This section will cover 1) general mechanism underlying PIR, 2) the effect of magnitude and duration of inhibitory current application, and 3) the effect of CaT and CaL currents on PIR.The default parameter values for this section are  +', = 20,  )*( = 1,  +'$ = 5,  *+& = 2.For the simulation of inhibitory input, we used the parameter  '--in the differential equation of membrane potential  (Eq.1).Default magnitude and duration of applied inhibitory current are -20 (that is  '--= −20) and 500ms.

Mechanisms underlying PIRs. Bifurcation analysis.
Figure 6A shows an example of PIR.Before a cell is given an inhibitory input, the cell exhibits a spontaneous tonic firing with a frequency of around 10 Hz.At t = 500 msec,  '-- was changed from 0 to -20 for 500ms and then the cell was hyperpolarized over this period.
Once the inhibitory input was removed at t = 1000 msec, the cell exhibited a burst where the frequency decreased over time while the magnitude increased.Although there is no clear way to define the duration of PIR, we can loosely define it as the time from the removal of inhibition to the moment when the trajectory returns to its original tonic firing solution (Fig. 6E-F).In terms of inter-spike intervals, this means that inter-spike interval returns to its original value in spontaneous tonic firing solution.In Fig. 6A, the trajectory returned to its original spontaneous tonic firing around t = 1600ms, thus the duration of PIR is around 600ms in this example.
In previous sections we observed that bursting solutions and spiking solutions are facilitated by the CaT current and modulated by the HCN-current when the CaL current is less significant.Therefore, we chose  and ℎ as bifurcation parameters and studied the underlying mechanism of PIRs. Figure 6B-D show bifurcation surfaces in (, ℎ, ) space with the projection of full model PIR solution shown in Fig. 6A.As before, green surface denotes -nullsurface for  +', = 20.Figure 6B shows the spontaneous tonic spiking solution until t = 500ms, which is similar to the one in Fig. 2B. Figure 6C illustrates what happens during the constant inhibitory input.Under inhibition, the surface of fixed points (S) moves rightward and as a result, the RK line of S is also shifted rightward.
Then, the trajectory jumps down to the lower part of S (lower white surface with black grid lines in the figure) and moves along the surface.Recall that  and ℎ increase at the same time along the surface.Let G denote the intersection curve of lower part of S and -nullsurface.This curve G is the set of globally stable fixed points and shown in the figure as a red line on the lower part of S. Since G lies on the lower surface of S, the trajectory (blue) moves along the lower surface of S to approach G and stays there until the inhibition is removed.In other words, this intersection curve G holds the trajectory until the removal of the inhibition and delimits the maximum level of  and ℎ during the inhibition.Hence, PIR duration is also delimited by this curve G. Once the inhibition is removed, the RK line goes back to its original place, the trajectory jumps up into the regime of stable periodic orbits, and a burst begins (Fig. 6D).While spiking within a burst, both  and ℎ decrease but ℎ decreases faster at the beginning of the burst.When the cell is released from inhibition, trajectory is close to the AH point so that the burst frequency is high.As the trajectory traverses the regime of stable periodic orbits toward the end of the regime (SNIC) on the left, the frequency of PIR decreases.After crossing the RK line, which is identical to the SNIC line (cf.Fig. 2), the trajectory approaches tonic spiking solution (Fig. 6B).
Figure 6E shows a two-parameter bifurcation diagram with the projection of the PIR solution (blue).The RK line is also shown as a black slant line.Thick blue curves near RK line correspond to spontaneous tonic spiking before the inhibition and after PIR.Once the inhibition is turned on,  and ℎ begin to increase (blue diagonal line).Figure 6F shows what happens during and after the inhibition.Figure 6E corresponds to the lower left corner of Fig. 6F.Once the inhibition is turned on, RK line is shifted right upward (black dotted line on the upper right corner; the red line denotes G,cf.Fig. 6C).Now the trajectory approaches G. Once the inhibition is removed, the RK line goes back to its original place (black solid line on the left side) and the trajectory jumps up into the regime of stable periodic orbits and a burst begins (cf.Fig. 6D).While spiking within a burst, the trajectory approaches the spontaneous tonic firing solution (cf.Fig. 6B,E).

Effect of magnitude and duration of inhibitory input on PIRs. Bifurcation analysis.
In this section, we studied the effect of magnitude and duration of inhibitory input on PIRs using a two-parameter bifurcation diagram.Figure 7A shows voltage profiles for progressively larger values of negative  '--( '--= -5, -10, -20, and -30 from top to bottom).The duration of inhibitory input is 500 ms for all four cases.As the magnitude of inhibitory input increases (input is more hyperpolarizing), the duration of PIR also increases.For more hyperpolarizing inputs such as  '--= -30, PIR shows an initially high frequency response and then the frequency of spiking slowly goes down to its base value.Figure 7B shows a two-parameter diagram with the projection of solutions for the same values of  '--as the time-series in Fig. 7A.The black curve at lower left corner is the RK line for  '--= 0. The remaining thin curves are intersection curves (Gs) described in previous section.These intersection curves (Gs) follow the same color code as the projected trajectories.As the magnitude of inhibitory input increases, G shifts to the upper right corner.When  '--= −5 or −10, input duration 500 ms was long enough so that trajectory reached G and stayed there.When  '--= −20 or -30, on the other hand, inhibitory input is effectively removed while the trajectory is still approaching G.Although there are two different scenarios, depending on magnitude and duration of inhibitory input, we see that, as the magnitude of inhibitory input increases, trajectory traverses farther away from the RK line (when  '--= 0) and the duration of PIR increases as a result.This fact also implies that, for large magnitude inhibitory input, a trajectory is closer to the AH point when it is released from inhibition.In other words, a large magnitude inhibitory input means a proximity to the AH point.Thus, we observe high frequency spiking at the beginning of PIR for large magnitude inhibitory input.
Figure 7C shows voltage profiles with four different input durations.Here  '--is fixed at -20.As the input duration increases, the duration of PIR increases.We also observe that two PIRs for input duration 500ms and 750ms look similar.Figure 7D shows a twoparameter diagram with the projection of solutions when input durations are the same as for the time-series in Fig. 7C.Termination (removal) of inhibitory input is denoted by dot in each trajectory.In all cases, trajectories approach the intersection curve G (thin magenta) over inhibition.When duration is 750ms, trajectory crosses G but it is trapped by that curve while ℎ increases.Hence trajectories move along G until inhibition is terminated.Figure 7D illustrates that PIR durations for input duration 500ms and 750 ms are similar because  values at release are similar.

Effect of CaT and CaL currents on PIRs. Bifurcation analysis.
Figure 8A shows the effect of  +', on PIRs using three  +', values, 15, 20, and 30 (from top to bottom) with fixed  +'$ = 5, input duration 500 ms, and  '--= −20.As shown in Fig. 1, the frequency of spontaneous tonic firing increases as  +', increases.Burst duration in PIR shows a slight but not significant increase as  +', increases.In the twoparameter diagram (Fig. 8B), the increase of  +', results in the horizontal stretch of the projected trajectories.The intersection curves (Gs) are also shifted to the right.This figure shows that maximum levels of ℎ values are similar in all cases.Since PIR duration can be estimated by the maximum level of ℎ when  *+& is fixed, this explains why there is only a minor change in PIRs as  +', increases.When  +', = 30, the trajectory is close to the AH point when it is released from inhibition.Due to this proximity to the AH point, PIR shows high frequency and small-magnitude spiking at the beginning of PIR (cf.Fig. 6).Now, we consider the effect of the CaL current on PIRs. Figure 8C shows PIRs for  +'$ = 5, 10, and 15 (from top to bottom) with fixed  +', = 20, input duration 500 ms, and  '--= −20.As  +'$ increases, we observe that the duration of PIR increases substantially.Note that the dynamics of the trajectory until jumping up into the regime of stable periodic orbits is facilitated by the CaT current and the HCN-current as discussed in Sections 3 and 4. In fact, when a cell is either in spontaneous tonic firing or hyperpolarized by inhibitory input, [] is almost constant at a very low level.Consequently, the activity patterns of these two states (either spontaneous tonic firing state before inhibitory input or silent state under hyperpolarizing input) can be explained by two slow variables  and ℎ as before.This is confirmed by Fig. 8D, which shows the projection of three solutions (for different  +'$ values) and the corresponding intersection curves in (, ℎ)-plane when cell is hyperpolarized.Both the projected trajectories and intersection curves for different values of  +'$ ( +'$ = 5, 10, 15) are not distinguishable.In other words, either the spontaneous tonic firing state before inhibitory input or silent state under hyperpolarizing input do not depend on  +'$ values when  +'$ is small.
However, when a cell is released from inhibition and jumps into the regime of stable periodic orbits, the CaL current begins to play an important role in PIR.We can explain the dynamics of trajectory over PIR using [], , and .Figure 8E shows the projection of trajectories in (, , []) -space for  +'$ = 5 (blue), 10 (red), and 15 (magenta).We omitted the RK surface and the AH surface, and plotted the SNPO/SNIC surface only in the middle for clarity (cf.Fig. 5D). Figure 8F shows the same trajectories with a SNPO/SNIC curve when [] = 0.05.Once a cell is released from the inhibitory input, [] initially increases ,and then overall level of [] decreases over spiking within a burst.The overall decrease of [] level is due to the decrease of spiking frequency over the time course of the active phase of burst.The values of  and  also decrease over spiking within a burst.For small  +'$ values, for example,  +'$ = 5 (red) and 10 (blue),  values are also relatively small (around 2.5 when  +'$ = 5 and around 5 when  +'$ = 10).In these two cases, active phase of burst is terminated when trajectory crosses the SNPO surface (Fig. 8F).Recall that in Fig. 5, we showed that the branches of stable periodic orbits lie above the -nullcline, hence the bursting solution moves leftward and jumps down to the lower branch of the S-shaped curve of fixed points .After jumping down, in this case, trajectory approaches spontaneous tonic firing solution (Fig. 8F).On the other hand, if  +'$ is sufficiently large, the bursting solution is not terminated by the SNPO because the SNPO has a negative  value over large  values (Fig. 8F).As shown in Fig. 3, since the lower end part of the branch of stable periodic orbits including the SNPO lies below -nullcline, averaged  ( bbbbb ) is between 0 and the RK.Thus, bursting solution moves leftward to approach  bbbbb .And while doing so, the [] level keeps decreasing and  bbbbb increases slowly (Fig. 3E). Figure 8F shows that the trajectory for  +'$ = 15 (magenta) turns around while spiking.Now the [] level becomes sufficiently low and the  value is sufficiently small.Then the PIR dynamics undergoes a transition to the dynamics of spontaneous tonic firing, which is driven by the CaT current and the HCN-current.

Discussion
In this study, we presented a conductance-based single-compartment model of an STN neuron, which plays an important role in the pathophysiology of the basal ganglia in Parkinson's disease.STN neurons exhibit characteristic activity patterns such as: a slow rhythmic firing [Beurrier et al., 1999;Bevan and Wilson, 1999], a calcium dependent postinhibitory rebound (PIR) bursts [Bevan et al., 2002a], and slow rhythmic bursting under sustained hyperpolarization [Beurrier et al., 1999].Recent experiments showed that interaction between the T-type calcium (CaT) current and the L-type calcium (CaL) current plays an important role in the generation of STN activity patterns [Beurrier et al., 1999;Bevan and Wilson, 1999;Bevan et al., 2002a].
The first single-compartment model of an STN neuron was developed by Terman et al.
[ Terman et al., 2002] and this model was able to generate PIR bursts with the CaT current.
Gillies and Willshaw [Gillies and Willshaw, 2006] developed a multi-compartment model, which contained the CaT current, the CaL current, and the HCN current and was able to generate characteristic activity patterns as stated above.However, the interaction of compartments in the model appears to be essential for its dynamical regimes.On the other hand, Hahn and McIntyre [Hahn and McIntyre, 2010] developed a single-compartment model that contained the CaT current and the CaL current but this model does not exhibit PIR burst nor slow rhythmic bursting under sustained hyperpolarization.
The STN model in this study is, to the best of our knowledge, the first single-compartment STN model that is able to generate characteristic activity patterns of STN neurons, especially activity patterns under hyperpolarization (hyperpolarization-induced bursts and PIR bursts).To investigate the roles and effects of these currents in rhythm generation we performed a bifurcation analysis using slow variables.We found that 1) the HCN current promotes single-spike activity patterns rather than bursting rhythms while nonetheless being an essential component for the bursting rhythms, 2) the CaT current enables STN cells to display characteristic firing patterns under hyperpolarization (hyperpolarizationinduced bursts and PIR bursts), and 3) the CaL current enriches and reinforces these bursting rhythms under hyperpolarization and PIR.

Roles of HCN, CaT, and CaL currents
Experimental results showed that the HCN current promotes single-spike activity patterns rather than bursting rhythms [Atherton et al., 2010].The bifurcation analysis of our model showed that the increase of maximal conductance  *+& (making the HCN current stronger) tends to yield a higher chance for spiking solution (Fig. 4G).This fact resulted from the proximity of a trajectory (projection of full model solution onto bifurcation diagram) to the RK line in the bifurcation diagram (cf.Fig. 2A,C, Fig. 4F,G).The RK line is the set of fold bifurcation points where the lower branch of stable fixed points turns around to become the middle branch of the unstable fixed points in the bifurcation diagram.Easier access to the RK line yields a higher chance for obtaining spiking solution in general.In a spontaneous tonic spiking solution, for example, the proximity of the trajectory to the RK line results in shorter inter-spike intervals and a higher frequency of spiking solution (Fig. 2C).Similarly, in a hyperpolarization-induced bursting solution this results in a higher frequency of bursting solution with fewer spikes within a burst (Fig. 4B,G).In conclusion, the larger availability of the HCN current renders means there is an easier access to the RK line, which facilitates a tonic spiking solution.
Our model showed that the CaT current is necessary for activity patterns under hyperpolarization (hyperpolarization-induced burst or post-inhibitory rebound (PIR) burst).This fact, in turn, indicates that the CaT current enables STN cells to generate various firing patterns under hyperpolarizing stimuli within the basal ganglia.The blocking or disrupting the CaT current may mute the emergence of rebound responses and hyperpolarization-induced rhythmic bursting solution.Our model also shows that the addition of the CaL current makes the response of an STN cell to inhibitory stimuli more prominent.In spontaneous tonic spiking solution there is an abrupt jump in frequency when the conductance of the CaL current becomes large enough.In hyperpolarizationinduced bursting rhythms or PIR bursts, the CaL current allows the cell to generate substantially longer bursting responses.To summarize, the synergistic interaction of the CaT current and the CaL current enables an STN cell to respond to hyperpolarizing stimuli in a salient way, and this fact may implicate the important roles of the CaT current and the CaL current in the pathophysiology of the basal ganglia in disorders, such as in Parkinson's disease, noted for elevated burstiness of STN neurons.

Bifurcation analysis and dynamical mechanisms of firing patterns.
Bifurcation analysis allowed us to understand the effect of three currents (CaT, CaL, HCN) considered in this study on activity patterns of an STN cell under specific conditions.The availability of a current affects the structure of bifurcation diagram of fast subsystem, and the dynamics of slow variables with respect to the resulting bifurcation diagram yields an explanation of the mechanism underlying a specific activity pattern.More specifically, we found that the generation of various activity patterns depends on several factors: the relative position of the S-shaped curve of fixed points with respect to stable periodic orbits and -nullcline, the place where the branches of stable periodic orbits terminate, and the existence of saddle-node bifurcation of periodic orbits (SNPO).
In this study we utilized four slow variables (ℎ, , , and []) for bifurcation analysis.Here, [] is the calcium concentration, ℎ =  *+& *  ,  =  +', *  ,  =   *  1 *  2 where  ,  ,  4,7 are gating variables and  *+& ,  +', ,  +'$ are maximal conductances for the HCN current, the CaT current, and the CaL current, respectively.In fact, these four slow variables are not sufficiently slow, so sometimes the projection of full model solution onto a bifurcation diagram shows some mismatch.In Figure 3D, for example, the projection of spiking solutions is not confined to the region inside two stable surfaces of stable orbits.If we make the four slow variables much slower, then we can resolve this mismatch, but the resulting activity patterns might not be physiologically realistic.Slow variables that were used in this study, however, were sufficiently slow so that we were able to obtain an insight into the underlying mechanism through bifurcation analysis.
Biophysicallly, the limit of very slow variables will lead to extremes in the firing patterns of a neuron (such as extremely long burst duration etc.).While one may argue that there is no sharp boundary between spiking and bursting activity, our results indicate that in the physiologically relevant regimes, the activity patterns exhibit specific bursting (which, if pushed to a limit in a mathematical consideration, leads to a coherent bifurcation structure).

Concluding remarks
As the only excitatory nucleus within the basal ganglia with strong pallidal and other inhibitory inputs, the fact that the STN is able to generate various bursting rhythms under hyperpolarization has an important implication in the pathophysiology of the basal ganglia.Synchronous beta oscillations within the basal ganglia is a hallmark of Parkinson's disease and has been associated with pathological symptoms related to movement [Brown, 2003;Hutchison, 2004;Kühn et al., 2004;Brown, 2007;Hammond et al., 2007;Mallet et al., 2008;Ray et al., 2008;Eusebio and Brown, 2009;Kühn et al., 2009;Park et al., 2010;Oswal et al., 2013;Stein and Bar-Gad, 2013;Ahn et al., 2015] Over the past decades, many theories have been developed with respect to the origin of the excessive beta rhythms within the basal ganglia (see Introduction).Two types of these theories focus on the role of the STN-GPe network and have been at the center of attention.In the first case, beta oscillations are generated in the cortex and the STN-GPe network in the basal ganglia has an ability to resonate or otherwise respond with oscillations at this frequency [see, for example, discussion in Stein and Bar-Gad, 2013].In the second case, the STN-GPe network by itself has an ability to generate beta oscillations and plays an important role in maintaining the beta rhythms independently or via thalamus to the cortex connection [Bevan et al., 2002b;Mallet et al., 2008;Merrison-Hort and Borisyuk, 2013].

Captions
Figure 1 Dependence of spontaneous tonic firing frequency on intrinsic parameters.We present the effects of (A) external constant input ( '--), (B) the CaT current ( +', ), (C) the HCN current ( *+& ), and (D) the CaL current ( +'$ ).In all four cases, the frequency increases as the current strength increases.In the panel (D), there is an abrupt increase in the frequency.Red plot in the same figure shows the result when the timescale of the dynamics of gating variable for the HCN current (f) was increased by 50%.and -30 (cyan).Black curve at lower left corner is the RK line for  '--= 0. Remaining thin curves are intersection curves of the lower surface of stable fixed points and nullsurface (Gs) introduced in the previous section for each value of  '--.Trajectory moves along the lower part of S to approach G under inhibition.These curves (Gs) use the same color code with the projection of full model solutions.(C) PIRs with four different input durations (100ms, 250ms, 500ms, and 750ms from top to bottom).(D) Two-parameter diagram with the projection of solutions (thick closed curves) when input durations are 100ms (blue), 250ms (red), 500ms (magenta), and 750ms (cyan).Black curve at lower left corner is the RK line and thin magenta curve at right upper corner is G. Termination of the inhibitory input is denoted by dot in each trajectory.

Figure
Figure 3A and 3B show bifurcation diagrams of the fast subsystem in (, )-space with the bifurcation parameter  for fixed  and [] values.Here, [] value is 0.05 in Fig. 3A and 0.18 in Fig. 3B.In each figure,  values are 2 (black), 6 (blue), and 10(red).Green curves are -nullcline when  +', = 20.There are several things to note on these bifurcation diagrams.First, we see that the lower and middle parts of the S-shaped curve of fixed points () remain almost the same over various  values.Especially the dependence of  value of the right knee (RK) of  ( F% ) on  values is negligible.This is clearly shown in Fig.3C, where the vertical lines denote the RK lines for [] = 0.05 (black) and [] = 0.18 (blue).Second, there are two different termination mechanisms of stable periodic orbits depending on the  values.As shown in the previous section, the family of stable periodic orbits emanates from  at the AH-point.For smaller  values such as 2, stable periodic orbits terminate in a saddle-node on an invariant circle (SNIC) bifurcation.On the other hand, for larger  values such as 6 and 10, the stable periodic orbits turn around at saddle-node bifurcation of periodic orbits (SNPO) to become unstable periodic orbits.The third thing to note is the relative position of -nullcline with respect to , especially with respect to the RK of .As seen in the figure, -nullcline lies above the RK for smaller [] values (Fig.3A).But -nullcline intersects at the middle and lower parts of  for larger [] values (Fig.3B).

Figure 3
Figure 3 (A -B) Bifurcation diagrams of fast subsystems in (, )-space with bifurcation parameter  when  and [] are fixed.Here, [] is fixed at 0.05 in Fig 3A and at 0.18 in Fig 3B.In each figure,  = 2 (black), 6 (blue), and 10 (red).Green curve is nullcline for  +', = 20.(C) Two-parameter bifurcation diagrams with the projection of full model spiking solutions.Two vertical lines are RK lines when [] = 0.05 (black) and 0.18 (blue).Dotted lines that emanate from RK lines are SNPO lines.Short horizontal lines denote the projection of full model solutions.For smaller  +'$ values, the projection of full model solutions lies on black RK lines (from bottom to top,  +'$ = 5, 10, 15).For larger  +'$ values, the full model solutions lie inside the regime of stable periodic orbits (from

Figure 6
Figure 6 (A) An example of post-inhibitory rebound (PIR) burst.At t = 500 msec,  '--was changed from 0 to -20 for 500ms.(B-D) Bifurcation surfaces in (, ℎ,  L )-space with the projection of full model PIR solution (blue) shown in (A).Green surface denotes -nullsurface for  +', = 20.(B) Spontaneous tonic spiking solution (from t = 0 ms to t = 500 ms).(C) The behavior of the trajectory under applied inhibitory current (from t = 500 ms to t = 1000 ms).The red curve on the lower surface of stable fixed points (S) is the intersection curve between this surface and -nullsurface (G).(D) Activity pattern once the cell is released from inhibition.The cell jumps into the regime of stable periodic orbits.(E) Two-parameter bifurcation diagram with the projection of the PIR solution (blue) near RK line (black).When inhibition is turned on, the trajectory was pushed away from the RK line.(F) Activity patterns during and after the inhibition.Once the inhibition is turned on, RK line is shifted right upward (black dotted line on the upper right corner).The red line denotes G. Plot (E) is a magnification of the lower left corner of plot (F).

Figure 7
Figure 7 The effect of magnitude and duration of inhibitory input on PIRs.(A) PIRs for  '--= -5, -10, -20, and -30 from top to bottom.(B) Two-parameter diagram with the

Figure 8
Figure 8 Effect of CaT and CaL currents on PIRs for input duration 500ms and  '--= −20.(A -B) Effect of the CaT current on PIRs.(A) PIRs for  +', = 15, 20, and 30 from top to bottom.(B) Two-parameter diagram with the projection of solutions (thick closed curves) when  +', = 15 (blue), 20 (red), and 30 (magenta).Black curve at lower left corner is the RK line for  '--= 0. Remaining thin curves are intersection curves of the lower surface of stable fixed points and -nullsurface (Gs) introduced in the previous section.Trajectory moves along the lower part of S to approach G under inhibition.These curves (Gs) use the same color code with the projection of full model solutions.(C -F) Effect of the CaL current on PIRs.(C) PIRs for  +'$ = 5, 10, and 15 from top to bottom.(D) Two-parameter diagram with the projection of solutions and the corresponding intersection curves over the inhibitory input. +'$ = 5 (blue), 10 (red), and 15 (magenta).(E) Projection of trajectories in (, , []) -space for the same values of  +'$ .SNPO/SNIC surface is also shown.(F) The same trajectories with SNPO/SNIC curve when [] = 0.05. Figure1

Table 1 .
Table 1 lists the resulting values of the kinetic parameters.Values of kinetic parameters.Units for each parameter values are shown in the first row except  <,A ,  <,!7 ,  <,A ,  <,!7 whose units are mM.
. Although whether the STN-GPe network generates the excessive beta rhythms in vivo in Parkinsonian brain or not is still uncertain, these two theories demonstrate the important role of the STN-GPe network in the excessive beta rhythmicity.In this context, the presently investigated dynamical mechanisms promoting the ability of an STN cell to generate bursting rhythms under either transient or sustained hyperpolarization may underlie excessively synchronous beta rhythms observed in Parkinsonian basal ganglia.Finally, we would like to note the growing interest in the adaptive Deep Brain Stimulation (DBS) of the STN in Parkinson's disease.The development of effective control of beta-band activity may benefit from the availability of a relatively simple STN model like the one considered here, which captures the major dynamical characteristics of STN activity.