# Analysis of the role of the low threshold currents I_{T} and I_{h} in intrinsic delta oscillations of thalamocortical neurons

^{1}Consejo Nacional de Investigaciones Científicas y Técnicas, Física Estadística e Interdisciplinaria, Centro Atómico Bariloche, San Carlos de Bariloche, Argentina^{2}Comisión Nacional de Energía Atómica and Consejo Nacional de Investigaciones Científicas y Técnicas, Centro Atómico Bariloche and Instituto Balseiro, San Carlos de Bariloche, Argentina

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 I_{T} 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 m_{T} and the recovering variable h_{T} of the calcium channel I_{T} is sufficient to generate low threshold oscillations in the delta band. We also explored the role of the hyperpolarization activated cationic current I_{h} in this reduced model and determine that, albeit not required, I_{h} 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

respectively; where *I _{Kleak}* =

*g*(

_{Kleak}*V*−

*E*)

_{Kleak}*S*,

*I*=

_{Naleak}*g*(

_{Naleak}*V*−

*E*)

_{Naleak}*S*are the leak currents,

*I*is the injected current,

_{inj}*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*and

_{T}*h*are the activation and inactivation variables respectively and

_{T}*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*) is a non-linear function of potential and calcium concentration;

_{o},Ca_{i}where *Ca _{o}* and

*Ca*are the extracellular and the intracellular concentrations of Ca

_{i}^{++}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 *g _{h}* is the maximum conductance,

*m*is the activation variable,

_{h}*E*is the reversal potential and

_{h}*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*) into a two-dimensional model (containing only

_{T}/dt*dV/dt*and

*dh*), we assumed an instantaneous activation of I

_{T}/dt_{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 Runge-Kutta integration algorithm with a 0.01 ms time step (Press et al., 1992).

**Figure 1. Dynamical mechanisms of oscillations in the I _{T}-Leaks model**.

**(A)**Frequency-current plots of oscillations induced by sustained injection of depolarizing current from a hyperpolarized stable potential (black) and injection of hyperpolarizing current from a depolarized stable potential (red) for the 3D I

_{T}-Leaks model (solid lines) and for the 2D I

_{T}-Leaks model (dashed lines). The inset shows magnification of the region of mismatch (hysteresis) between the two plots at the lowest frequencies for the 3D model.

**(B)**Bifurcation diagram of the I

_{T}-Leaks model using the default value of

*p*(7.0 × 10

_{T}^{−5}cm/s) while maintaining other parameters at default. The diagram shows voltage at fixed points (black V/I curve) and max/min of limit cycle as the current injection is changed (dots). The dashed region on the V/I curve (black) and the red dots represent instability whereas the continuous line and blue dots represent regions of stable equilibrium and stable periodic orbits respectively.

**(C)**Bifurcation diagram after shifting the activation variable of I

_{T}by −3 mV, using

*p*= 3.0 × 10

_{T}^{−5}cm/s.

**(D)**Bifurcation diagram of the I

_{T}-leaks model with kinetic parameters and voltage dependence as in McCormick and Huguenard (1992).

*p*set to 11 × 10

_{T}^{−4}cm/s.

**(E)**Bifurcation diagram of the I

_{T}-leaks model using a

*p*value of 9.0 × 10

_{T}^{−5}cm/s while maintain other parameters at default.

**(F)**Bifurcation diagram after shifting m

_{T}by −3 mV, using a

*p*value of 4.0 × 10

_{T}^{−5}cm/s. Conventions for panels

**(C–F)**as in

**(B)**.

## Results

### 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*), the propensity to fire repetitive bursts is strongly affected by the permeability of the I

_{Naleak}_{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*= 5 × 10

_{T}^{−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*value of I

_{T}_{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*at 5 × 10

_{T}^{−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*, or a global depolarizing shift larger than +2 mV of the inactivation variable

_{T}*h*) favor sustained oscillations (see Figure 8 in Amarillo et al., 2014). Similarly, increasing the maximum permeability

_{T}*p*by less than 30% (from 5 to 7 × 10

_{T}^{−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*. 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

_{T}*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*set to 1.1 × 10

_{T}^{−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*to slightly higher values (above 8.0 × 10

_{T}^{−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*= 4.0 × 10

_{T}^{−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*the model should undergo a

_{T}*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*to 9.0 × 10

_{T}^{−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 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).

**Figure 2. Phase-plane portrait of the I _{T}-leaks model after reduction of dimensionality from 3D to 2D**.

**(A)**Using a

*p*value of 7.0 × 10

_{T}^{−5}cm/s, at a depolarized equilibrium (

*I*= 6 pA) the membrane potential is stable at −61.5 mV (blue trace in inset). At this potential, the

_{inj}*V*nullcline (blue line) intersects the

*h*nullcline (gray line) on a single point (

_{T}*a*). Decreasing

*I*to 2 pA induces oscillations (red trace in inset) and in the phase-plane plot, shifts the

_{inj}*V*nullcline upwards (red line) with loss of stability of the intersecting point (

*b*) and appearance of a stable limit cycle (supercritical Hopf bifurcation). Black orbits represent the variation of

*V*as function of

*h*in the direction indicated by the arrow heads.

_{T}**(B)**Using the same

*p*value as in

_{T}**(A)**, at a hyperpolarized equilibrium (

*I*= −7 pA) the membrane potential stabilized at −75.2 mV (blue trace in inset), the

_{inj}*V*(blue line) and

*h*(gray line) nullclines intersect in a single stable point (

_{T}*a*). Changing

*I*to −6 pA induces oscillations (red trace in inset) and, in the phase plane portrait, shifts the

_{inj}*V*nullcline (red line) downwards with loss of stability of the intersecting point (

*b*) and disappearance of an unstable limit cycle (subcritical Hopf bifurcation).

**(C)**After increasing

*p*to 9.0 × 10

_{T}^{−5}cm/s, at a hyperpolarized equilibrium (

*I*= −11 pA), the membrane potential is stabilized at −77.7 mV (blue trace in inset). In the phase plane portrait,

_{inj}*V*(blue line) and

*h*(gray line) nullclines intersect in three points; one stable point (

_{T}*a*) and two unstable points (

*b*and

*c*). At a lower level of hyperpolarization (

*I*= −10 pA), depolarizing current shifts the

_{inj}*V*nullcline (red line) downwards with destruction of the points

*a*and

*b*, leaving one unstable point (

*d*) and allowing the model to oscillate at very low frequencies. Bottom panel shows a magnification corresponding to the encircled area on the upper panel. Conventions for panels

**(B,C)**as in

**(A)**.

### 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*). Injection of hyperpolarizing current—up to −6 pA—further removes the inactivation of I

_{T}_{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.

**Figure 3. I _{h} has different effects on the two bifurcations**.

**(A)**Upper trace shows the time course of the membrane potential (V) during oscillations elicited at −1 pA (near the upper bifurcation) in the absence of I

_{h}. The time courses of I

_{T}and its gating variables are shown aligned underneath. The bottom diagram shows the time course of the relative contribution of the three currents considered (outward currents in blue and inward currents represented by shades of red and yellow). The white solid line separates total outward from total inward components and therefore represents the net depolarizing or hyperpolarizing drive of the model at a given time point. Vertical dotted lines are positioned at: the peak of the oscillation (a); the valley of the oscillation (b); and the time point of maximum contribution of I

_{T}(c). Horizontal dotted lines represent zero values for currents and gating variables and 50% of the total current in diagrams of relative contribution.

**(B)**Time course of I

_{T}, its gating variables and the relative contribution of the currents during oscillations elicited at −4.5 pA (near the lower bifurcation) in the absence of I

_{h}.

**(C)**Time course of I

_{T}, its gating variables and the relative contribution of the currents for a large hyperpolarizing current injection (without oscillations) in the absence of I

_{h}. Vertical dotted line is at the onset of current injection.

**(D)**Time course of I

_{T}, I

_{h}, their gating variables and the relative contribution of the currents during oscillations elicited at −3 pA (near the upper bifurcation) after switching on I

_{h}.

**(E)**Time course of I

_{T}, I

_{h}, their gating variables and the relative contribution of the currents during oscillations produced under strong hyperpolarization (−23 pA, near the lower bifurcation) in the presence of I

_{h}. Conventions in panels

**(B–E)**as in panel

**(A)**.

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*has to be increased above 1.5 × 10

_{T}^{−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).

**Figure 4. I _{h} adds robustness to the oscillations**. Bifurcation diagram of the I

_{T}-I

_{h}-Leaks model using a

*p*value of 7 × 10

_{T}^{−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*= −8 pA (upper trace) and oscillations elicited in the presence of Ih and absence of spiking mechanisms using the same value of

_{inj}*I*(lower trace). The bars above the bifurcation diagram indicate the range of

_{inj}*I*that elicits oscillations in the model with spiking mechanisms both in the presence of I

_{inj}_{h}(blue bar) and in the absence of I

_{h}(gray bar).

## 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 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 and implemented 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*for limit cycle—gray dots in Figure 4 with the range of

_{inj}*I*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

_{inj}^{+}-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*that elicit oscillations does not change (from −2 to −31 pA; compare the range of

_{inj}*I*for limit cycle—blue dots in Figure 4 with the range of

_{inj}*I*that elicit oscillations when spiking mechanisms are present—blue bar above the bifurcation diagram on Figure 4).

_{inj}### 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 structure—and 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 saddle-node 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 saddle-node 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.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We thank Dr. Inés Samengo for helpful discussions. This work was subsidized by grant PIP 112 200901 00738 (CONICET-Argentina).

## References

Amarillo, Y., and Nadal, M. S. (2013). “Dynamic analysis of the minimal requirements for intrinsic rhythmic burst firing in thalamocortical neurons,” in *Program. 615.629/F618, Society for Neuroscience Meeting 2013* (San Diego, CA).

Amarillo, Y., Zagha, E., Mato, G., Rudy, B., and Nadal, M. S. (2014). The interplay of seven subthreshold conductances controls the resting membrane potential and the oscillatory behavior of thalamocortical neurons. *J. Neurophysiol*. 112, 393–410. doi: 10.1152/jn.00647.2013

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Anderson, M. P., Mochizuki, T., Xie, J., Fischler, W., Manger, J. P., Talley, E. M., et al. (2005). Thalamic Cav3.1 T-type Ca2+ channel plays a crucial role in stabilizing sleep. *Proc. Natl. Acad. Sci. U.S.A*. 102, 1743–1748. doi: 10.1073/pnas.0409644102

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Budde, T., Caputi, L., Kanyshkova, T., Staak, R., Abrahamczik, C., Munsch, T., et al. (2005). Impaired regulation of thalamic pacemaker channels through an imbalance of subunit expression in absence epilepsy. *J. Neurosci*. 25, 9871–9882. doi: 10.1523/JNEUROSCI.2590-05.2005

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Crunelli, V., and Hughes, S. W. (2010). The slow (<1 Hz) rhythm of non-REM sleep: a dialogue between three cardinal oscillators. *Nat. Neurosci*. 13, 9–17. doi: 10.1038/nn.2445

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Destexhe, A., Neubig, M., Ulrich, D., and Huguenard, J. (1998). Dendritic low-threshold calcium currents in thalamic relay cells. *J. Neurosci*. 18, 3574–3588.

Destexhe, A., and Sejnowski, T. J. (2003). Interactions between membrane conductances underlying thalamocortical slow-wave oscillations. *Physiol. Rev*. 83, 1401–1453. doi: 10.1152/physrev.00012.2003

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ermentrout, B. (2002). *Simulating, Analyzing, and Animating Dynamical Systems: A Guide to Xppaut for Researchers and Students*. Philadelphia, PA: SIAM. doi: 10.1137/1.9780898718195

Ernst, W. L., Zhang, Y., Yoo, J. W., Ernst, S. J., and Noebels, J. L. (2009). Genetic enhancement of thalamocortical network activity by elevating alpha 1g-mediated low-voltage-activated calcium current induces pure absence epilepsy. *J. Neurosci*. 29, 1615–1625. doi: 10.1523/JNEUROSCI.2081-08.2009

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hansel, D., Mato, G., and Meunier, C. (1995). Synchrony in excitatory neural networks. *Neural Comput* 7, 307–333. doi: 10.1162/neco.1995.7.2.307

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hansel, D., Mato, G., and Pfeuty, B. (2012). “The role of intrinsic cell properties in synchrony of neurons interacting via electrical synapses,” in *Phase Response Curves in Neuroscience*, eds N. W. Schultheiss, A. A. Prinz, and R. J. Butera (New York, NY: Springer), 361–398. doi: 10.1007/978-1-4614-0739-3_15

Hughes, S. W., Cope, D. W., and Crunelli, V. (1998). Dynamic clamp study of Ih modulation of burst firing and delta oscillations in thalamocortical neurons *in vitro*. *Neuroscience* 87, 541–550. doi: 10.1016/S0306-4522(98)00170-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hutcheon, B., Miura, R. M., Yarom, Y., and Puil, E. (1994). Low-threshold calcium current and resonance in thalamic neurons: a model of frequency preference. *J. Neurophysiol*. 71, 583–594.

Izhikevich, E. (2005). *Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*. Cambridge, MA: MIT Press.

Jahnsen, H., and Llinas, R. (1984). Ionic basis for the electro-responsiveness and oscillatory properties of guinea-pig thalamic neurones *in vitro*. *J. Physiol*. 349, 227–247. doi: 10.1113/jphysiol.1984.sp015154

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kim, D., Song, I., Keum, S., Lee, T., Jeong, M. J., Kim, S. S., et al. (2001). Lack of the burst firing of thalamocortical relay neurons and resistance to absence seizures in mice lacking alpha(1G) T-type Ca(2+) channels. *Neuron* 31, 35–45. doi: 10.1016/S0896-6273(01)00343-9

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lee, J., Kim, D., and Shin, H. S. (2004). Lack of delta waves and sleep disturbances during non-rapid eye movement sleep in mice lacking alpha1G-subunit of T-type calcium channels. *Proc. Natl. Acad. Sci. U.S.A*. 101, 18195–18199. doi: 10.1073/pnas.0408089101

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lee, S.-C., Cruikshank, S. J., and Connors, B. W. (2010). Electrical and chemical synapses between relay neurons in developing thalamus. *J. Physiol*. 588, 2403–2415. doi: 10.1113/jphysiol.2010.187096

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Leresche, N., Lightowler, S., Soltesz, I., Jassik-Gerschenfeld, D., and Crunelli, V. (1991). Low-frequency oscillatory activities intrinsic to rat and cat thalamocortical cells. *J. Physiol*. 441, 155–174. doi: 10.1113/jphysiol.1991.sp018744

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lewis, T. J., and Skinner, F. K. (2012). “Understanding activity in electrically coupled networks using PRCs and the theory of weakly coupled oscillators,” in *Phase Response Curves in Neuroscience*, eds N. W. Schultheiss, A. A. Prinz, and R. J. Butera (New York, NY: Springer), 329–359. doi: 10.1007/978-1-4614-0739-3_14

Ludwig, A., Budde, T., Stieber, J., Moosmang, S., Wahl, C., Holthoff, K., et al. (2003). Absence epilepsy and sinus dysrhythmia in mice lacking the pacemaker channel HCN2. *EMBO J*. 22, 216–224. doi: 10.1093/emboj/cdg032

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lytton, W. W., Destexhe, A., and Sejnowski, T. J. (1996). Control of slow oscillations in the thalamocortical neuron: a computer model. *Neuroscience* 70, 673–684. doi: 10.1016/S0306-4522(96)83006-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mancilla, J. G., Lewis, T. J., Pinto, D. J., Rinzel, J., and Connors, B. W. (2007). Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. *J. Neurosci*. 27, 2058–2073. doi: 10.1523/JNEUROSCI.2715-06.2007

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mato, G., and Samengo, I. (2008). Type I and type II neuron models are selectively driven by differential stimulus features. *Neural Comput*. 20, 2418–2440. doi: 10.1162/neco.2008.10-07-632

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

McCormick, D. A., and Bal, T. (1997). Sleep and arousal: thalamocortical mechanisms. *Annu. Rev. Neurosci*. 20, 185–215. doi: 10.1146/annurev.neuro.20.1.185

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

McCormick, D. A., and Feeser, H. R. (1990). Functional implications of burst firing and single spike activity in lateral geniculate relay neurons. *Neuroscience* 39, 103–111. doi: 10.1016/0306-4522(90)90225-S

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

McCormick, D. A., and Huguenard, J. R. (1992). A model of the electrophysiological properties of thalamocortical relay neurons. *J. Neurophysiol*. 68, 1384–1400.

McCormick, D. A., and Pape, H. C. (1990). Properties of a hyperpolarization-activated cation current and its role in rhythmic oscillation in thalamic relay neurones. *J. Physiol*. 431, 291–318. doi: 10.1113/jphysiol.1990.sp018331

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ortuno, T., Grieve, K. L., Cao, R., Cudeiro, J., and Rivadulla, C. (2014). Bursting thalamic responses in awake monkey contribute to visual detection and are modulated by corticofugal feedback. *Front. Behav. Neurosci*. 8:198. doi: 10.3389/fnbeh.2014.00198

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Pfeuty, B., Mato, G., Golomb, D., and Hansel, D. (2003). Electrical synapses and synchrony: the role of intrinsic currents. *J. Neuroscience* 23, 6280–6294.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1992). *Numerical Recipes in FORTRAN 77*, Vol. 1. New York, NY: Press Syndicate of the University of Cambridge.

Rush, M. E., and Rinzel, J. (1994). Analysis of bursting in a thalamic neuron model. *Biol. Cybern*. 71, 281–291. doi: 10.1007/BF00239616

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Samengo, I., Mato, G., Elijah, D. H., Schreiber, S., and Montemurro, M. A. (2013). Linking dynamical and functional properties of intrinsically bursting neurons. *J. Comp. Neurosci*. 35, 213–230. doi: 10.1007/s10827-013-0449-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Santoro, B., Chen, S., Lüthi, A., Pavlidis, P., Shumyatsky, G. P., Tibbs, G. R., et al. (2000). Molecular and functional heterogeneity of hyperpolarization-activated pacemaker channels in the mouse CNS. *J. Neurosci*. 20, 5264–5275.

Soltesz, I., Lightowler, S., Leresche, N., Jassik-Gerschenfeld, D., Pollard, C. E., and Crunelli, V. (1991). Two inward currents and the transformation of low-frequency oscillations of rat and cat thalamocortical cells. *J. Physiol*. 441, 175–197. doi: 10.1113/jphysiol.1991.sp018745

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Traub, R. D., Buhl, E. H., Gloveli, T., and Whittington, M. A. (2003). Fast rhythmic bursting can be induced in layer 2/3 cortical neurons by enhancing persistent Na+ conductance or by blocking BK channels. *J. Neurophysiol*. 89, 909–921. doi: 10.1152/jn.00573.2002

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Keywords: T-type calcium channel, thalamocortical neurons, repetitive burst firing, sub-threshold conductances

Citation: Amarillo Y, Mato G and Nadal MS (2015) Analysis of the role of the low threshold currents I_{T} and I_{h} in intrinsic delta oscillations of thalamocortical neurons. *Front. Comput. Neurosci*. **9**:52. doi: 10.3389/fncom.2015.00052

Received: 12 February 2015; Accepted: 21 April 2015;

Published: 07 May 2015.

Edited by:

Andre Longtin, University of Ottawa, CanadaReviewed by:

Frances K. Skinner, University of Toronto, CanadaSteven James Cox, Rice University, USA

Copyright © 2015 Amarillo, Mato and Nadal. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Yimy Amarillo, Consejo Nacional de Investigaciones Científicas y Técnicas, Física Estadística e Interdisciplinaria, Centro Atómico Bariloche, Avenida Bustillo 9500, San Carlos de Bariloche, R8402AGP Rio Negro, Argentina, amarillo@cab.cnea.gov.ar