Analysis of the role of the low threshold currents IT and Ih in intrinsic delta oscillations of thalamocortical neurons

Thalamocortical neurons are involved in the generation and maintenance of brain rhythms associated with global functional states. The repetitive burst firing of TC neurons at delta frequencies (1–4 Hz) has been linked to the oscillations recorded during deep sleep and during episodes of absence seizures. To get insight into the biophysical properties that are the basis for intrinsic delta oscillations in these neurons, we performed a bifurcation analysis of a minimal conductance-based thalamocortical neuron model including only the IT channel and the sodium and potassium leak channels. This analysis unveils the dynamics of repetitive burst firing of TC neurons, and describes how the interplay between the amplifying variable mT and the recovering variable hT of the calcium channel IT is sufficient to generate low threshold oscillations in the delta band. We also explored the role of the hyperpolarization activated cationic current Ih in this reduced model and determine that, albeit not required, Ih amplifies and stabilizes the oscillation.


Introduction
Repetitive burst firing of thalamocortical (TC) neurons in the delta band has been linked to the expression of the rhythms that characterize slow wave sleep and the pathological spike and wave discharges of absence epilepsy (McCormick and Bal, 1997;Budde et al., 2005). The synchronized expression of repetitive bursting in TC neurons of the behaving animal is the result of the interaction between the intrinsic properties of these neurons and the synaptic activity of the thalamo-reticulo-cortical network (Lytton et al., 1996;Destexhe and Sejnowski, 2003). The prominent role of intrinsic ionic mechanisms in the generation and maintenance of the oscillations at the cellular level has been extensively demonstrated: individual TC neurons maintained in vitro are able to fire bursts repetitively either spontaneously or after injection of hyperpolarizing current, and this ability is indeed conserved under conditions of synaptic isolation (McCormick and Pape, 1990;Leresche et al., 1991). In addition, the EEG expression of rhythms associated with repetitive burst firing of TC neurons (slow wave oscillations during deep sleep and spike and wave discharges during absence seizures) are strongly affected by genetic or pharmacological manipulation of the ion channels expressed by TC neurons (Kim et al., 2001;Ludwig et al., 2003;Lee et al., 2004;Anderson et al., 2005;Budde et al., 2005). Specifically, genetic elimination of the calcium channel pore forming subunit CaV3.1 (the main channel subunit carrying I T in TC neurons from mice) abolishes the generation of low threshold spikes (LTS) (Kim et al., 2001), and suppresses the delta oscillations during NREM sleep (Lee et al., 2004). Conversely, overexpression of this channel subunit results in a phenotype of pure absence epilepsy in mice (Ernst et al., 2009). An absence epilepsy phenotype is also obtained by genetic elimination of HCN2 (Ludwig et al., 2003), which is an I h channel subunit strongly expressed in TC neurons (Santoro et al., 2000).
It has been previously shown that the two firing regimes of thalamocortical neurons, tonic and burst firing, can be described by two distinct oscillatory systems that operate independently at different membrane potentials and at different time scales (Rush and Rinzel, 1994). In a previous study, we characterized the seven different conductance's that operate at hyperpolarized membrane potentials, and established their contributions to intrinsic delta oscillations. In that study, we showed that the minimal ionic mechanisms required for generation and maintenance of oscillations compatible with physiological (or pathological) repetitive burst firing are I T and the leak currents (Amarillo et al., 2014). Here we analyze the bifurcation structure of this minimal model and show that the system can enter the oscillatory regime (limit cycle) from two different stable equilibriums (which occur at physiologically plausible membrane potentials) via different dynamic mechanisms: the transition from a depolarized stable equilibrium occurs via a supercritical Hopf bifurcation, whereas at hyperpolarized potentials the system can undergo either a subcritical Hopf bifurcation or a saddle-node bifurcation on invariant cycle (Izhikevich, 2005). We discuss possible functional, physiological and pathophysiological implications of this dynamic behavior.
Although a role of I h in repetitive burst firing of TC neurons has been suggested previously (McCormick and Pape, 1990;Soltesz et al., 1991;Hughes et al., 1998), our simulations indicate that I h is not essential for repetitive burst firing (Amarillo et al., 2014). It has been previously demonstrated that the main role of I h is the stabilization of the resting membrane potential (RMP), and that this current is one of the main determinants of the positive shift of the RMP from the potassium equilibrium potential in TC neurons (Amarillo et al., 2014 and references therein). In this study, we use phase plane and bifurcation analysis to determine the role of Ih in repetitive burst firing and we show that I h stabilizes the oscillations by increasing the voltage range (and the range of current injection magnitudes) at which stable oscillations occur. Some of these results have been presented previously in abstract form (Amarillo and Nadal, 2013).

Methods
The HH-like equations used in this study have been published elsewhere (Amarillo et al., 2014). Briefly, the voltage equations for the minimal I T -Leaks model and the I T -I h -Leaks model are: and  Table 2.  Table 2). respectively; where I Kleak = g Kleak (V−E Kleak )S, I Naleak = g Naleak (V−E Naleak )S are the leak currents, I inj is the injected current, C is the membrane capacitance and S the surface of the neuron (see parameter values in Tables 1, 2).
The equations used for I T are: where p T is the maximum permeability, m T and h T are the activation and inactivation variables respectively and m T∞ (V), h T∞ (V), τ mT (V), τ hT (V) are the steady state and time constants of activation and inactivation (see Tables 1, 2). G(V, Ca o ,Ca i ) is a non-linear function of potential and calcium concentration; where Ca o and Ca i are the extracellular and the intracellular concentrations of Ca ++ and z, F, R, and T are the valence, the Faraday constant, the gas constant and the absolute temperature respectively. The equations for I h are: whereḡ h is the maximum conductance, m h is the activation variable, E h is the reversal potential and m h∞ (V) and τ mh (V) are the steady state activation and time constant respectively (see Table 1). Bifurcation analysis and phase plane portrait analysis were performed using XPPAUT (Ermentrout, 2002). In order to convert the three-dimensional I T -leaks model containing three differential equations (dV/dt, dm T /dt, and dh T /dt) into a twodimensional model (containing only dV/dt and dh T /dt), we assumed an instantaneous activation of I T and replaced the time dependent equation for the steady state equation of the m variable of I T . Thus, the current Equation (3) was replaced by Frequency current plots ( Figure 1A) were obtained by averaging the inter LTS intervals occurring during each of 4000 current steps of 10 s, between −10 and +10 pA, using a second order Condition Global voltage shifts are applied on all voltage dependence parameters of a given gating variable of I T as previously performed (Destexhe et al., 1998;Amarillo et al., 2014). *In Figures 1E,F the parameters of voltage dependence are as indicated in the table, however the values of p T are 9 × 10 −5 and 4 × 10 −5 cm/s respectively (see text). Frontiers in Computational Neuroscience | www.frontiersin.org Runge-Kutta integration algorithm with a 0.01 ms time step (Press et al., 1992).

Phase Plane and Bifurcation Analysis of a Minimal Model of Thalamocortical Neurons
In a previous study, we used a murine TC neuron conductance model to explore the minimal requirements for generating periodic oscillations at delta frequencies (Amarillo et al., 2014). In the minimal model, which only has the calcium current I T and the leak currents (I Kleak and I Naleak ), the propensity to fire repetitive bursts is strongly affected by the permeability of the I T current. The baseline value of maximum permeability of I T in this minimal model corresponds to the experimental value obtained from rodent TC neurons in vitro: p T = 5 × 10 −5 cm/s. In agreement with the rodent experimental data, the TC model has a low propensity to fire repetitive bursts when using this baseline p T value of I T . However, increasing the availability of I T -by manipulating its maximum permeability or its voltage dependence of activation and inactivation-enables periodic low threshold oscillations. Thus, while keeping the value of p T at 5 × 10 −5 cm/s, changes in the gating variables of I T that result in a larger window current component (a global hyperpolarizing shift larger than −2 mV in the activation variable m T , or a global depolarizing shift larger than +2 mV of the inactivation variable h T ) favor sustained oscillations (see Figure 8 in Amarillo et al., 2014). Similarly, increasing the maximum permeability p T by less than 30% (from 5 to 7 × 10 −5 cm/s) induces spontaneous oscillations of 32 mV of amplitude at 2.3 Hz. We take this value (7 × 10 −5 cm/s) as default in the rest of this paper.
Here we analyze the transitions from two stable equilibriums, one occurring at a relatively depolarized membrane potential V (positive to about −63 mV) and the other occurring at hyperpolarized V (negative to about −75 mV) into a stable limit cycle by using injected current I inj as parameter. We first analyzed the bifurcation structure of the 3-dimentional system comprising the differential Equations (1), (4), and (5). Figure 1A shows frequency current plots of the system using the default value of p T . These plots reveal a small range of bistability (hysteresis, see inset) which is bounded in the left side by a fold limit cycle bifurcation and on the right side by a subcritical Hopf bifurcation. In this bifurcation, the stable limit cycle coalesces with an unstable limit cycle and both disappear. We also compared the 3D dynamical system with a 2D system obtained by making instantaneous the kinetics of the m gate of I T (see Methods), i.e., making m T = m T∞ (V). This is justified because there is a one order of magnitude difference between activation and inactivation time constants (with activation being faster than inactivation). We found that the change of dimensionality from 3D to 2D does not modify the oscillatory activity of the model in response to either injection of hyperpolarizing current from a positive equilibrium or injection of depolarizing current from a negative equilibrium. In fact we found that the bifurcation structure is exactly the same, although the positions of the bifurcations points are changed ( Figure 1A).
The bifurcation diagram shows that the transition from the resting state to the limit cycle and from the limit cycle to rest occurs via a subcritical Hopf bifurcation at hyperpolarized potentials ( Figure 1B, open arrow) whereas the transitions at depolarized voltages occur via a supercritical Hopf bifurcation ( Figure 1B, filled arrow).
We repeated this analysis using different voltage dependences (global shifts) of the activation and inactivation gates of I T and found the same bifurcation structure. Figure 1C shows the bifurcation diagram of the I T -Leaks model after shifting the activation of I T by −3 mV using p T = 3.0 × 10 −5 cm/s. Furthermore, the bifurcation diagram of the I T -Leaks model using the same values of voltage dependence of activation of I T as in the seminal study by McCormick and Huguenard (1992) (see Table 2), with p T set to 1.1 × 10 −4 cm/s, also have the same structure ( Figure 1D). This indicates that the dynamical properties of I T are insensitive to small variations of voltage dependence provided that the maximum permeability of I T is kept to the minimum required to enable sustained oscillations. However, setting p T to slightly higher values (above 8.0 × 10 −5 cm/s) changes the type of bifurcation occurring at hyperpolarized potentials from a subcritical Hopf to a saddle-node bifurcation on invariant cycle, as the I/V relationship changes from monotonic to non-monotonic ( Figure 1E). Similarly, setting p T = 4.0 × 10 −5 cm/s after shifting the activation of I T by -3 mV, also changes the type of bifurcation occurring at hyperpolarized potentials from subcritical Hopf to saddle node on invariant cycle ( Figure 1F). This dynamical behavior indicates that at some intermediate value of p T the model should undergo a Bogdanov-Takens codimension-2 bifurcation (Izhikevich, 2005). On the other hand the bifurcation occurring at depolarized potentials stays a supercritical Hopf for all values of p T .
To explore this further, we analyzed the 2D model using a graphic visualization in phase plane portraits. We analyzed the phase plane portraits near the bifurcation points using either the default value of p T (Figures 2A,B) or a slightly increased value (9.0 × 10 −5 cm/s, Figure 2C). Figure 2A shows the transition from a depolarized equilibrium to the limit cycle as hyperpolarizing current is injected beyond the threshold for generating repetitive LTSs. At equilibrium, the nullclines intersect at a single stable point (labeled a). Injection of hyperpolarizing current displaces the V nullcline upwards and destabilizes the intersecting point through a supercritical Hopf bifurcation initiating the oscillation (superimposed orbits). Figure 2B shows the transition from a hyperpolarized equilibrium to the limit cycle after injection of depolarizing current injection. At equilibrium, the nullclines intersect each other at a single stable point. Injection of depolarizing current beyond the LTSs threshold shifts the V nullcline downwards; this destabilizes the intersecting point through a subcritical Hopf bifurcation, thereby initiating the oscillation. The same transition is shown in Figure 2C after increasing p T to 9.0 × 10 −5 cm/s. In this case, the nullclines intersect at three points at equilibrium. Injection of depolarizing current beyond the LTS threshold shifts the V nullcline downwards and cause the destruction of two of the intersecting Frontiers in Computational Neuroscience | www.frontiersin.org points (points a and b), triggering the abrupt appearance of the oscillation. The change of the dynamical behavior from one type of bifurcation to another-that in this case depends both on the magnitude of current injection and on the maximum permeability of I T -indicates again that the model can undergo a codimension-2 Bogdanov-Takens bifurcation, which is the only possible codimension-2 bifurcation in 2-dimensional systems (Izhikevich, 2005).

Understanding the Role of I H in Periodic Burst Firing
We next explored the interaction between I h and I T in the presence of the leak conductances. With p T set to the default value and I h switched off, oscillations occur in a narrow range of current injection values (approximately between +2 and −6 pA; Figure 1B). With no current or minimal values of current injection (Figure 3A), the availability of I T is small because the inactivation is large (small h T ). Injection of hyperpolarizing current-up to −6 pA-further removes the inactivation of I T and gives rise to oscillations with larger amplitudes (Figure 3B). Under these conditions, oscillations are maintained solely by the regenerative activation of I T . Injection of hyperpolarizing current larger than −6 pA induces stronger de-inactivation of I T ; yet, activation never develops because the drive of the current injection is overpowering ( Figure 3C). With I h switched on (Figures 3D,E), the RMP is shifted toward depolarized values due to the steady activation of I h . In this case, injection of low values of hyperpolarizing current produces a similar sequence of de-inactivation/activation of I T as in the absence of I h (compare I T gating variables in Figures 3A,D). Notice that both the current magnitude and the gating variables of I T reach similar values with or without little activation of I h . In contrast, larger magnitudes of hyperpolarizing current (negative to −6 pA) lead to stronger activation of I h , which results in the removal of a larger fraction of inactivation of I T (Figure 3E). Hence, the depolarizing drive contributed by I h not only adds to the regenerative activation of I T during the ascending phase of the LTS, giving rise to larger and faster oscillations, but it also permits the occurrence of these oscillations at much hyperpolarized levels (in contrast to the case shown in Figure 3C). The resulting effect is that I h strongly affects the bifurcation that occurs at hyperpolarized potential yet it affects weakly the bifurcation that takes place at depolarized potential.
The bifurcation diagram of this four-dimensional system [with differential Equations (2), (4), (5) and (8)] shows an I/V relationship that is monotonic for a large range of p T values. In order to change the I/V relationship to non-monotonic under these conditions, p T has to be increased above 1.5 × 10 −4 cm/s. This does not imply, however, that a saddle-node bifurcation takes place because for the values of the current where the I/V curve becomes non-monotonic, the resting state has already lost stability via a Hopf bifurcation (data not shown). The diagram shows that the system enters the limit cycle via a supercritical Hopf bifurcation at depolarized potentials, and via a subcritical Hopf at hyperpolarized potentials (Figure 4). The bifurcation diagram also shows that the range of current injection that elicits oscillations is much broader in the presence of I h than in the minimal I T -leaks model (Figure 4), consistent with a role of I h in stabilizing the oscillations. This effect of Ih is independent of the sodium spiking mechanism (Figure 4, blue and gray bar; see Discussion).

Discussion
Thalamocortical neurons have two regimes of excitability: (1) the tonic firing mode that occurs at membrane potentials positive to about −60 mV and characterized by firing of action potentials at frequencies that correlates linearly with the strength of the stimuli (McCormick and Feeser, 1990); and (2) the firing of stereotyped bursts of action potentials at high frequency at membrane potentials negative to about −65 mV (Jahnsen and Llinas, 1984). Tonic firing is said to be a relay mode that allows information to be reliably transmitted to the targeted cortical areas (i.e., TC neurons under tonic mode function as integrators) and consistently occurs during active cognitive states. The role of burst firing is less clear; it has been proposed that TC neuron bursting could play a role in FIGURE 4 | I h adds robustness to the oscillations. Bifurcation diagram of the I T -I h -Leaks model using a p T value of 7 × 10 −5 cm/s while maintaining other parameters at default values (see Table 1 and Amarillo et al., 2014). Superimposed in gray is the max/min of limit cycle of the I T -leaks model (see Figure 1B). Vertical dotted lines labeled (A-E) are positioned at the values of current injection used in the corresponding panels of Figure 3. Conventions are as in Figure 1. The inset shows an example of the repetitive bursting oscillations elicited in the presence of I h and the spiking mechanisms (I Na and I K ) using I inj = −8 pA (upper trace) and oscillations elicited in the presence of Ih and absence of spiking mechanisms using the same value of I inj (lower trace). The bars above the bifurcation diagram indicate the range of I inj that elicits oscillations in the model with spiking mechanisms both in the presence of I h (blue bar) and in the absence of I h (gray bar). stimulus detectability and/or in switching the cortical targets from inattentive states to the activated states that characterize focused attention (Weyand et al., 2001;Ortuno et al., 2014). On the other hand, periodic bursting of TC neurons correlates with brain states characterized behaviorally by cognitive arrest (deep sleep, absence seizures), and physiologically, by global synchronization of the thalamocortical system (i.e., TC neurons under repetitive burst firing mode function as resonators).
In an early study by Rush and Rinzel (1994), the bimodal excitability of TC neurons was explained by the co-existence of two oscillatory systems that operate independently at different membrane potentials and at different temporal scales. A fast system, composed by the fast amplifying sodium conductance and the fast resonant potassium conductances, operates at depolarized potentials and underlies tonic firing, whereas a slow system, composed by the gating variables (fast amplifying activation and slow resonant inactivation) of I T , operates at hyperpolarized membrane potentials and underlies LTS. At the most depolarized phase of these LTSs, the fast system is activated generating the characteristic burst of action potentials. It has been shown, both experimentally with TTX (see Figure 10 in McCormick and Pape, 1990) and computationally (see Figure 2C in Rush and Rinzel, 1994), that the oscillatory behavior of TC neurons is unaffected by the sodium spikes (see inset in Figure 4). These evidences support the idea that the two oscillatory systems can be studied separately and that the fast system does not affect the behavior of the slow system. Indeed, we compared the response to current injection of the minimal I T -leaks model with and without spiking mechanisms (HH-like models of fast Na+ and K+ currents taken from Traub et al., 2003 andimplemented as in Amarillo et al., 2014), and found no differences in the range of I inj that elicit oscillations (from −6 to +2 pA; compare the range of I inj for limit cycle-gray dots in Figure 4 with the range of I inj that elicit oscillations when spiking mechanisms are present-gray bar above the bifurcation diagram in Figure 4). The only difference between the two sets of simulations is the presence of fast bursts of Na + -K + spikes riding on the LTSs that reached the spike threshold when spiking mechanisms are present. Furthermore, we made a similar comparison with and without spiking mechanisms in the I T -I h -leaks model, and the range of I inj that elicit oscillations does not change (from −2 to −31 pA; compare the range of I inj for limit cycle-blue dots in Figure 4 with the range of I inj that elicit oscillations when spiking mechanisms are present-blue bar above the bifurcation diagram on Figure 4).

Neurocomputational Properties of the Minimal I T -Leaks Model
In the present study, we focus on the dynamic structure of the slow system (given by the gating variables of I T ) that underlies repetitive burst firing of TC neurons. The bifurcation structureand the upper bifurcation in particular-of the minimal model that supports this slow system is consistent with the structure of a resonator, as it has been previously suggested (Hutcheon et al., 1994). The supercritical Hopf bifurcation occurring as the system enters a limit cycle from depolarized potentials enables oscillations of graded amplitude. Thus, minor hyperpolarizations induce oscillations of low amplitude without the requirement of a threshold potential (in contrast to the full blown oscillations characteristic of integrators, which usually require a large displacement of membrane potential to reach a threshold). These low amplitude oscillations could be potentiated by synaptic inputs arriving at the same frequency (resonant frequency), eventually reaching the threshold of the fast oscillatory system. The consequent action potential firing leads to neurotransmitter release, and therefore, to inter-neuronal communication, thereby promoting the synchronization of the thalamocortical system. We propose that TC neurons could anomalously enter the repetitive burst firing mode via this supercritical Hopf bifurcation during activated (depolarized) states. This would result in the pathological synchronization of the thalamocortical system at delta frequencies during wakefulness, which correlates with the occurrence of episodes of absence epilepsy.
Our analysis also shows that the behavior of the transition from a stable equilibrium to the limit cycle at hyperpolarized membrane potentials is much more complex. Depending on the value of the maximal I T permeability, the system can undergo either a subcritical Hopf bifurcation or a saddlenode on invariant cycle bifurcation. This means that the neurocomputational behavior of the system can change from a resonator to an integrator, and, significantly, that this switch is controlled by the level of maximal I T permeability. In addition, the time scale of the resonator increases continuously as the Bogdanov-Takens co-dimension 2 bifurcation is approached (Izhikevich, 2005).
This flexibility could have some interesting consequences on the functionality of TC cells, since resonators and integrators are driven by different optimal stimuli. Integrators fire in response to scale-free depolarizing stimuli (i.e., stimuli whose time scale depends on the firing rate of the neuron but does not depend on times scales of the intrinsic dynamics). Resonators, instead, are most efficiently driven by input stimuli containing both depolarizing and hyperpolarizing phases, with significant power in the frequency band corresponding to the intrinsic frequencies of the cell (Mato and Samengo, 2008). Choosing the adequate value of the I T permeability would allow the system to control its sensitivity to inputs in the delta band or in the infra-slow band (Crunelli and Hughes, 2010) (by approaching the Bogdanov-Takens bifurcation). These properties could be very important in situations where bursting is not perfectly periodic. In Samengo et al. (2013) it was shown that, given the adequate level of variability, bursters are able to codify input information and that the coding mechanism is essentially determined by the bifurcation structure.
The other important consequence of the transition from resonator to integrator lies in the collective network behavior. Both types of neuronal behavior tend to have very different synchronization properties. Resonators tend to synchronize when these neurons interact with chemical excitatory interactions that have a time constant shorter than the period of oscillation; whereas the same interaction has a desynchronizing effect on networks of integrators (Hansel et al., 1995). For inhibition, the relation is the inverse, at least for not very strong values of the coupling constant.
I T permeability can also affect network dynamics via its effect on gap junctions. This type of intercellular communication has been found in thalamocortical cells in early postnatal stages (Lee et al., 2010). The effect of gap junctions can be modulated by intrinsic currents (Pfeuty et al., 2003;Hansel et al., 2012). For instance, Pfeuty et al. (2003) showed that changing the values of sodium and potassium conductances allows to control the degree of network synchrony, from fully synchronized to a completely asynchronous behavior. In Mancilla et al. (2007) it is shown that this effect can be the opposite in some cases (see Lewis and Skinner, 2012 for a discussion on the discrepancy). In any case, the modulation of I T permeability could also affect developmental processes via neuronal synchronization.
In summary, there are several mechanisms that would permit to control the dynamical state of the network and the flow of information through the thalamocortical system just by regulating up or down the permeability of I T .
The Role of I H I h not only contributes to the stabilization of the RMP in TC neurons (Amarillo et al., 2014), but it also adds robustness to the low threshold oscillations by potentiating the initial phase of depolarization and also by allowing larger excursions of the membrane potential between LTSs (favoring de-inactivation of I T ). Oscillations elicited at hyperpolarized potentials require the activation of I h to provide a recovering depolarization (pacemaker potential), which adds to the depolarizing influence of a residually activated I T . Loss of the stabilizing effects of I h alters both the robustness of the oscillations and the voltage regime at which they occur, giving rise to more readily, yet aberrant, oscillations at more depolarized potentials. This conclusion is also supported by the absence epilepsy phenotype of the I h KO mice (HCN2 principal subunit; Ludwig et al., 2003), which would otherwise contradict an essential role of I h in repetitive burst firing.
In addition to these effects on the rhythmic bursting behavior of TC neurons, introduction of I h tends to suppress the saddlenode bifurcation that is present in the minimal I T -leaks model. Hence, the presence of I h would favor the resonator behavior. This means that a more complex modulation would be required to switch the system to the integrator behavior: i.e., concomitant down-regulation of I h and up-regulation of I T .
The present study contributes to unveil the complex dynamic behavior of TC neurons. Our conclusions are in line with the suggestion that these cells have the potential to perform a sophisticated role in controlling and processing the information that flows from sensory sources to the cortex (primary thalamic nuclei), and between different cortical areas via higher order thalamic nuclei.