Analysis of the Neuron Dynamics in Thalamic Reticular Nucleus by a Reduced Model

Strategically located between the thalamus and the cortex, the inhibitory thalamic reticular nucleus (TRN) is a hub to regulate selective attention during wakefulness and control the thalamic and cortical oscillations during sleep. A salient feature of TRN neurons contributing to these functions is their characteristic firing patterns, ranging in a continuum from tonic spiking to bursting spiking. However, the dynamical mechanism under these firing behaviors is not well understood. In this study, by applying a reduction method to a full conductance-based neuron model, we construct a reduced three-variable model to investigate the dynamics of TRN neurons. We show that the reduced model can effectively reproduce the spiking patterns of TRN neurons as observed in vivo and in vitro experiments, and meanwhile allow us to perform bifurcation analysis of the spiking dynamics. Specifically, we demonstrate that the rebound bursting of a TRN neuron is a type of “fold/homo-clinic” bifurcation, and the tonic spiking is the fold cycle bifurcation. Further one-parameter bifurcation analysis reveals that the transition between these discharge patterns can be controlled by the external current. We expect that this reduced neuron model will help us to further study the complicated dynamics and functions of the TRN network.

1. Introduction

The thalamic reticular nucleus (TRN) is a brain area containing a large population of GABAergic neurons (Sherman and Guillery, 2001; Pinault, 2004; Halassa and Acsády, 2016; Crabtree, 2018). It receives glutamatergic projections from the thalamocortical and corticothalamic neurons (Sherman and Guillery, 2001; Pinault, 2004; Crabtree, 2018), and inhibitory and modulatory inputs from many subcortical regions (Jourdain et al., 1989; Govindaiah et al., 2010; Beierlein, 2014; Herrera et al., 2016; Nakajima et al., 2019). However, all its efferents are targeted on cells in the thalamic nuclei (Sherman and Guillery, 2001; Pinault, 2004; Crabtree, 2018), and nearly every part of the thalamus receives inhibition from the TRN (Swanson et al., 2019). Located in such a strategical position, TRN has long been hypothesized to play an important role in the regulation of information flow transferred through the thalamus. In an influential theoretical proposal, Francis Crick suggested that “if the thalamus is the gateway to the cortex, the reticular complex might be described as the guardian of the gateway” (Crick, 1984). This hypothesis of the importance of TRN on selective attention has been confirmed by recent studies on rodents, monkeys, and humans (McAlonan et al., 2006, 2008; Halassa et al., 2014; Ahrens et al., 2015; Wimmer et al., 2015; Wells et al., 2016; Nakajima et al., 2019). Moreover, a growing body of experiments starts to draw the conclusion that TRN is also a center for the control of the brain oscillation during sleep. The pioneering work of Steriade et al. (1985, 1987) found that, in cats in vivo, the TRN-disconnected thalamic nuclei lose their ability to generate spindle oscillations, while the isolated TRN deprived of inputs from the cortex and thalamus can itself generate spindle rhythms, demonstrating TRN is the spindle pacemaker. Moreover, optogenetic studies have shown that activation pulses on TRN neurons can evoke spindle oscillations in the connected thalamic or cortical areas (Halassa et al., 2011; Barthó et al., 2014; Clemente-Perez et al., 2017). Another study in awake mice found that local activation of TRN can rapidly induce delta oscillation in the corresponding region of the cortex area, along with a decrease in arousal state (Lewis et al., 2015). Putting together, these experimental results suggest that TRN has strong functional controls on attentional filtering and oscillation regulation in the brain.

A key ingredient of TRN neurons to exert their functions in the brain is their rich firing patterns, including bursts (low-threshold Ca²⁺ potentials) when neurons are hyperpolarized (e.g., during sleep) (Steriade et al., 1986; Huguenard and Prince, 1992; Bal and McCormick, 1993; Llinás and Steriade, 2006), tonic spikes when neurons are depolarized (e.g., during the wakefulness or attention state) (Steriade et al., 1986, 1993; Herd et al., 2013; Rovó et al., 2014; Lewis et al., 2015), and firing patterns between them (Pinault, 2004; Halassa and Acsády, 2016). On the one hand, the intrinsic bursting in TRN neurons is important for the regulation of brain oscillations. Barthó et al. (2014) demonstrated that the number of spikes in each TRN burst predicts the progression of a spindle event. Moreover, genetic deletion of Ca²⁺ channels in TRN neurons in mice abolished low-threshold Ca²⁺ potentials and bursts, which further suppressed spindle rhythms while enhanced delta oscillations during sleep (Pellegrini et al., 2016; Fernandez et al., 2018; Vantomme et al., 2019; Li et al., 2020). On the other hand, TRN tonic spikes are prevalent during the wakefulness or attentional state (Pinault, 2004; McAlonan et al., 2006, 2008). Such tonic inhibition permits information transfer from the thalamus during an attentional focus (McAlonan et al., 2006, 2008; Coulon and Landisman, 2017), and allows fine-grained inhibitory control on the thalamic targets (McAlonan et al., 2008; Wimmer et al., 2015; Halassa and Acsády, 2016).

To reproduce these characteristic firing patterns, detailed conductance-based neuron models were proposed for TRN neurons (Destexhe et al., 1996; Bazhenov et al., 1998, 2002; Vijayan and Kopell, 2012; Fan et al., 2017). However, in these full models, too many voltage-dependent conductance channels and dynamical parameters are involved, which makes it hard to analyze the neuronal dynamics clearly and, hence, prevents us from understanding the inner mechanism underlying the rich dynamics of TRN neurons. Furthermore, it has been documented that distributed TRN cells in different regions usually exhibit different spiking properties (Brunton and Charpak, 1997; Lee et al., 2007; Clemente-Perez et al., 2017; Li et al., 2020; Martinez-Garcia et al., 2020). A full neuron model with too many parameters makes it hard to identify the key elements that contribute to this intrinsic heterogeneity in firing patterns. Without understanding the dynamics of single TRN neurons clearly, it will be difficult for us to understand the much more complicated dynamics and functions of the TRN network.

In this study, motivated by the foregoing shortcomings of a full neuron model, we present a reduced three-variable model for TRN neurons by adopting a reduction method proposed by Kepler et al. (1992). The main idea of this reduction algorithm is to group different variables operating in a similar time scale into a single representative variable. The reduced model is an approximation of the original full model but preserves its key dynamical features. The reduced model also keeps the original biological meanings of the corresponding full model, as the individual channel currents can be recovered from the reduced variables. The key advantage of the reduced TRN neuron model is that we can perform theoretical analysis on its dynamics. Utilizing the reduced model, we carry out an analysis to elucidate the bifurcation mechanisms under two characteristic discharge patterns of TRN neurons, i.e., the rebound bursting and tonic spiking. We show that the rebound bursting is a type of “fold/homo-clinic” bifurcation and the tonic spiking is the fold cycle bifurcation. We also carry out a one-parameter bifurcation analysis and find that transition between these two firing patterns can be achieved by varying the external input. We hope that this study will pave the way for us to study the much more complicated dynamics of the TRN network in the future study.

2. A Full Model of TRN Neurons

The full mathematical model of TRN neurons we consider was proposed by Bazhenov et al. (1998, 2002), which is described as a single-compartment Hodgkin-Huxley schema, with a set of active channels sufficient to produce the typical firing patterns seen in TRN neurons. Similarly to other dynamical systems that describe isopotential excitable neural membranes, this model could be written as this form:

$\begin{array}{ll}C\frac{\text{d}V}{\text{d}t}+I(V,\left\{x_i\right\})=\frac{10^{-3}{I}_{syn}(t)}{A},& (1)\end{array}$

where V is the membrane potential, C = 1 μF/cm² is the membrane capacitance, I(V, {x_i}) is a term expressed as a function of V and the gating variables x_i describing the ion currents, and I_syn is the received synaptic current which is normalized by the membrane area A = 1.43 × 10⁻⁴cm² (Bazhenov et al., 2002). Specifically, this model contains three active ionic currents, which are fast sodium current I_Na and a delayed rectifier I_K for spiking generation, and a low-threshold Ca²⁺ current I_T for rebound bursting, and two leaky currents, which are a membrane leakage current I_L and a potassium leaky current I_KL controlled by neuromodulators like acetylcholine and norepinephrine (McCormick, 1992). The ion current term is written as,

$\begin{array}{ll}I(V,\left\{x_i\right\})=I_L+I_{KL}+I_{Na}+I_K+I_T,& (2)\end{array}$

where the leakage current is modeled as I_L = g_L(V−E_L), with the leakage channel conductance $g_L=0.06\text{ mS}/\text{c}{\text{m}}^2$ and the leakage reversal potential E_L = −70 mV, and the potassium leaky current is modeled as I_KL = g_KL(V − E_KL), with the reversal potential E_KL = −100 mV. The last three items are voltage-dependent ionic currents which are expressed in a unified form as,

$\begin{array}{ll}{I}_i=g_i{M}^j{N}^k(V-E_i),i\in \left\{Na,K,T\right\}& (3)\end{array}$

where g_i is the maximal conductance, M is the activation gating variable, N is the inactivation gating variable, j and k are the numbers of the gating variables, and E_i is the reversal potential. Specifically, $I_{Na}=g_{Na}{m}^{3}h(V-E_{Na})$, $I_K=g_K{n}^4(V-E_K)$, and $I_T=g_T{p}^{2}q(V-E_T)$, with the maximal conductance $g_{Na}=100\text{ mS}/\text{c}{\text{m}}^2$, $g_K=10\text{ mS}/\text{c}{\text{m}}^2$, and the reversal potentials E_Na = 50 mV, E_K = −100 mV, and E_T = 120 mV. The dynamics of gating variables are expressed in the below unified form,

$\begin{array}{ll}\frac{d\theta }{dt}={\varphi }_{\theta }\frac{{\theta }_{\infty }(V)-\theta }{{\tau }_{\theta }(V)},\theta \in \left\{m,h,n,p,q\right\},& (4)\end{array}$

where θ_∞ is the voltage-dependent steady state, τ_θ is the voltage-dependent time constant, and ϕ_θ is a constant. For a gating variable θ ∈ {m, h, n}, θ_∞ and τ_θ are expressed by the voltage-dependent state transition variables α and β, i.e.,

$\begin{array}{ll}{\theta }_{\infty }(V)=\frac{{\alpha }_{\theta }(V)}{{\alpha }_{\theta }(V)+{\beta }_{\theta }(V)},{\tau }_{\theta }=\frac{1}{{\alpha }_{\theta }(V)+{\beta }_{\theta }(V)}.& (5)\end{array}$

The equations and values of the above variables are listed in Table 1.

TABLE 1

Table 1. Kinetics of gating variables.

3. Reduction of the Full Model

As described above, the full model has six dynamical variables (V, m, h, n, p, and q), making it hard to analyze the inner mechanism and difficult to accommodate TRN neurons with different firing patterns. In this study, by adopting a reduction method proposed for the conductance-based neuron model (Kepler et al., 1992), we present a reduced model with a minimal set of variables to capture the behaviors of the original full model. The overall reduction procedure can be summarized into the following four steps.

3.1. Step 1: Converting Gating Variables Into Equivalent Potentials

In the original expressions, variables are not dimensionally equivalent. To combine different variables for reduction, making them dimensionally equivalent is necessary. Membrane potential provides a connection between these variables, so we consider that all gating variables are converted into equivalent potentials v_θ, in term of that in the voltage-clamp recording, equilibrium voltage will gives rise to the value θ, i.e., θ = θ_∞(v_θ) or equivalently, $v_{\theta }={\theta }_{\infty }^{-1}(\theta )$.

Using the equivalent potentials, the original full model Equation (1) can be re-expressed as,

$\begin{array}{ll}C\frac{\text{d}V}{\text{d}t}+F(V,\left\{v_{\theta }\right\})=\frac{10^{-3}{I}_{syn}(t)}{A},& (6)\end{array}$

where

$\begin{array}{l}F(V,\left\{v_{\theta }\right\})=I(V,\left\{{\theta }_{\infty }(v_{\theta })\right\})\\ =g_{Na}{m}_{\infty }^3(v_m)h_{\infty }(v_h)\\ +g_T{p}_{\infty }^2(v_p)q_{\infty }(v_q)(V-E_T)\\ +g_K{n}_{\infty }^4(v_n)(V-E_K)+I_L+I_{KL}.& (7)\end{array}$

From Equation (5), we also obtain

$\begin{array}{ll}\frac{d{v}_{\theta }}{dt}={\varphi }_{\theta }\frac{{\theta }_{\infty }(V)-{\theta }_{\infty }(v_{\theta })}{{\tau }_{\theta }(V){\theta }_{\infty }^{\prime }(v_{\theta })}=f_{\theta }(V,v_{\theta }),\theta \in \left\{m,h,n,p,q\right\},& (8)\end{array}$

where ${\theta }_{\infty }^{\prime }$ denotes the derivative of θ_∞ with respect to v_θ.

3.2. Step 2: Grouping Variables at a Similar Time Scale

By simulation, we obtain the dynamics of equivalent potentials in the full TRN neuron model, which are displayed in Figure 1. We see that they can be categorized into three classes. The first category includes V and v_m, which have the largest variations. The second category includes v_h, v_n, and v_p, which have an intermediate amplitude of variation and change on a slower time scale compared to the first category. Note that v_h and v_n are anti-synergistic to v_p, in terms that the increases in v_h and v_n hyper-polarize the membrane potential while the increase in v_p depolarizes the membrane potential. The third category only includes v_q, which has the smallest amplitude of variation and changes at the slowest speed. These properties suggest that the dynamics of a TRN neuron in the six-dimensional space can be effectively reduced to a three-dimensional manifold.

FIGURE 1

Figure 1. The evolution of the membrane potential and the equivalent potentials in the full thalamic reticular nucleus (TRN) neuron model. The model (Equations 6–8) receives an inhibitory current with the size of −0.03 mS/cm² during the first 200 ms. In the next 400 ms, the model exhibits a rebound bursting. V and v_m have the fastest time scale and exhibit the largest variations; v_h, v_n, and v_p response at a slower time scale and display intermediate amplitude variations; v_q has the slowest time constant and shows the modest variation. The neuronal parameters used are E_L = –77 mV, $g_{KL}=0.00793\text{ mS}/\text{c}{\text{m}}^2$, and $g_T=2\text{ mS}/\text{c}{\text{m}}^2$.

To locate this manifold, we group those equivalent potentials behaving similarly into a single representative variable as a weighted average of their values, which are written as,

$\begin{array}{ll}x={\rho }_{V}V+{\rho }_m{v}_m,& (9)\end{array}$

$\begin{array}{ll}y={\rho }_h{v}_h+{\rho }_n{v}_n+{\rho }_p{v}_p,& (10)\end{array}$

$\begin{array}{ll}z=v_q,& (11)\end{array}$

where all the coefficients {ρ_i}, for i ∈ {V, m, h, n, p, q}, are non-negative, and they satisfy ρ_V + ρ_m = 1 and ρ_h + ρ_n + ρ_p = 1.

3.3. Step 3: Determining the New Variables

We fix the coefficients for the new variables in Equations (9–11) under the condition that the dynamics of the membrane potential of the reduced model can approximate that of the full model as accurately as possible.

From Equation (9), we have

$\begin{array}{ll}\frac{dx}{dt}={\rho }_V\frac{dV}{dt}+{\rho }_m\frac{d{v}_m}{dt}.& (12)\end{array}$

Substituting Equations (6–8) into Equation (12) and approximating the dynamics in the first-order Taylor expansion, we have,

$\begin{array}{l}C\frac{\text{d}x}{\text{d}t}={\rho }_V\frac{{10}^{-3}{I}_{syn}(t)}{A}\\ -[{\rho }_{V}F(x,y,z)+{\delta }_h{\rho }_V\frac{\partial F}{\partial v_h}+{\delta }_n{\rho }_V\frac{\partial F}{\partial v_n}+{\delta }_p{\rho }_V\frac{\partial F}{\partial v_p}\\ +{\delta }_V\left({\rho }_V\frac{\partial F}{\partial V}-\frac{C{\varphi }_m{\rho }_m}{{\tau }_m}\right)\\ +{\delta }_m\left({\rho }_V\frac{\partial F}{\partial v_m}+\frac{C{\varphi }_m{\rho }_m}{{\tau }_m}\right)+O\left({\delta }^2\right)],& (13)\end{array}$

where δ_j = x−j, for j ∈ {V, v_m}, and δ_l = y−l, for l ∈ {v_h, v_n, v_p}. It is straightforward to check that ρ_hδ_vh + ρ_nδ_vn + ρ_pδ_vp = 0, and ρ_Vδ_V + ρ_mδ_vm = 0. To get the above result, we have used the conditions that f_m(x, x) = 0 and ∂f_m(x, x)/∂V = −∂f_m(x, x)/∂v_m = ϕ_θ/τ_m.

We determine the coefficients {ρ_i} by requiring that the discrepancy between the reduced and the full models (i.e., the first-order Taylor expansion in Equation 13) be as small as possible. This is formulated as an optimization problem. Specifically, for determining the coefficients ρ_h, ρ_n, ρ_p, ρ_m, andρ_V, the optimization problem is formulated as follows,

$\begin{array}{l}min{\displaystyle |{\delta }_{v_h}{\rho }_h\frac{\partial F}{\partial v_h}+{\delta }_{v_n}{\rho }_n\frac{\partial F}{\partial v_n}+{\delta }_{v_p}{\rho }_{v_p}\frac{\partial F}{\partial v_p}}\\ +{\delta }_V\left({\rho }_V\frac{\partial F}{\partial V}-\frac{C{\varphi }_m{\rho }_m}{{\tau }_m}\right)\\ {\displaystyle +{\delta }_{v_m}\left({\rho }_V\frac{\partial F}{\partial v_m}+\frac{C{\varphi }_m{\rho }_m}{{\tau }_m}\right)|}\\ \begin{array}{cccc}\text{s}.\text{t}.& {\rho }_h{\delta }_{v_h}+{\rho }_n{\delta }_{v_n}+{\rho }_p{\delta }_{v_p}=0,& & {\rho }_h,{\rho }_n,{\rho }_p\ge 0,\\ & {\rho }_V{\delta }_V+{\rho }_m{\delta }_{v_m}=0,& & {\rho }_m,{\rho }_V\ge 0,\end{array}& (14)\end{array}$

where the symbol |·| denotes the absolute value. This optimization problem can be solved by using the Lagrangian method (as shown in details in Supplementary Material), which gives,

$\begin{array}{ll}{\rho }_p=k,& (15)\end{array}$

$\begin{array}{ll}{\rho }_h=(1-k)\frac{\partial F}{\partial v_h}/\left(\frac{\partial F}{\partial v_h}+\frac{\partial F}{\partial v_n}\right),& (16)\end{array}$

$\begin{array}{ll}{\rho }_n=(1-k)\frac{\partial F}{\partial v_n}/\left(\frac{\partial F}{\partial v_h}+\frac{\partial F}{\partial v_n}\right),& (17)\end{array}$

$\begin{array}{ll}{\rho }_V=\frac{-(C\frac{{\varphi }_m}{{\tau }_m}+\frac{\partial F}{\partial V})+\sqrt{{(C\frac{{\varphi }_m}{{\tau }_m}+\frac{\partial F}{\partial V})}^2-4C\frac{{\varphi }_m}{{\tau }_m}(\frac{\partial F}{\partial V}+\frac{\partial F}{\partial v_m})}}{-2(\frac{\partial F}{\partial V}+\frac{\partial F}{\partial v_m})},& (18)\end{array}$

$\begin{array}{ll}{\rho }_m=1-{\rho }_V,& (19)\end{array}$

where 0 < k < 1 is an undetermined constant, which we normally set k to be 0.01.

3.4. Step 4: Formulating the Reduced Model

After fixing the new variables, we can write down their approximated dynamics. Specifically, for the variable x (corresponding to the membrane potential), its dynamics is written as,

$\begin{array}{l}\frac{C}{{\rho }_V}\frac{dx}{dt}=-g_{Na}{m}_{\infty }^3(x)h_{\infty }(y)(x-E_{Na})-g_K{n}_{\infty }^4(y)(x-E_K)\\ -g_T{p}_{\infty }^2(y)q_{\infty }(z)(x-E_T)-I_L-I_{KL}+\frac{10^{-3}}{A}{I}_{syn}(t).& (20)\end{array}$

In this study, we have ignored the discrepancy term in Equation (13) from the full model. Note that x has become an alternative variable of membrane potential V (in the following section we refer to V as x).

For the variable y, based on Equation (10), we get its dynamics,

$\begin{array}{l}\frac{dy}{dt}={\rho }_h\frac{d{v}_h}{dt}+{\rho }_n\frac{d{v}_n}{dt}+{\rho }_p\frac{d{v}_p}{dt},\\ ={\rho }_h{\varphi }_h\frac{h_{\infty }(x)-h_{\infty }(y)}{{\tau }_h(x)h_{\infty }^{\prime }(y)}+{\rho }_n{\varphi }_n\frac{n_{\infty }(x)-n_{\infty }(y)}{{\tau }_n(x)n_{\infty }^{\prime }(y)}\\ +{\rho }_p{\varphi }_p\frac{p_{\infty }(x)-p_{\infty }(y)}{{\tau }_p(x)p_{\infty }^{\prime }(y)}.& (21)\end{array}$

In this study, we have used the approximations v_h ≈ y, v_n ≈ y, and v_p ≈ y.

For the variable z, based on Equation (11), its dynamics is written as,

$\begin{array}{ll}\frac{dz}{dt}={\varphi }_q\frac{q_{\infty }(x)-q_{\infty }(z)}{{\tau }_q(x)q_{\infty }^{\prime }(z)}.& (22)\end{array}$

The above Equations (20–22) form the reduced model of TRN neurons. The details of the model are given in Supplementary Material: The details of the reduced model.

4. Analysis of the Reduced Model

In the above, we have reduced the original six-variable TRN neuron model into a dynamical system with three variables. In this section, we investigate the firing patterns and perform a bifurcation analysis of the reduced model. All simulations and analyses are performed using the BrainPy library (Wang et al., 2021).

4.1. Firing Patterns in the Reduced Model

We first verify that our reduced TRN neuron model behaves similarly to the original full one and can reproduce the firing patterns as seen in experiments. Two models are set to have the same parameters, including the neuronal parameters and the external input currents (Figure 2). Under a constant excitatory input with the size of 0.2 μA/cm², two models first burst, then gradually produce the same tonic spikes (Figure 2A). After shortening the depolarization duration (the excitatory current is only applied during the first 200 ms, Figure 2B), the reduced model also shows the same oscillation behavior. Specifically, both models first burst, then display tonic spikes, and finally exhibit subthreshold oscillations. Moreover, injection of a hyperpolarizing current pulse (an inhibitory current with the size of −0.05 μA/cm² and the duration of 200 ms) results in a similar rebound bursting (Figure 2C). Therefore, the reduced model can effectively capture the TRN neuron firing patterns, including tonic and burst firings, as seen in TRN neurons in vivo and in vitro (Contreras et al., 1992; Bal and McCormick, 1993; Lee et al., 2007).

FIGURE 2

Figure 2. Firing pattern comparison between the original (Equation 1) and the reduced (Equations 20–22) TRN neuron model. (A) Two models are stimulated with a constant excitatory input 0.2 μA/cm². They both exhibit decelerated burst and stationary tonic spikes. (B) Two models receive a step current with the size of 0.2 μA/cm². The removal of the external current results in a similar subthreshold oscillation in both models. (C) An inhibitory current with the size of −0.05 μA/cm² and the duration of 200 ms is applied to two models. Both rebound bursting and subthreshold oscillation observed in the full model are reproduced in the reduced one. The parameters used in this simulation are g_T = 2.25 mS/cm², g_KL = 0.0152 mS/cm², and k = 0.01.

In the next, we are going to apply dynamical system theory to illustrate the bifurcation mechanism underlying such firing patterns. We consider a fast-slow decomposition (Rinzel, 1985; Rinzel and Lee, 1987) based on the time scale differences between three variables. Specifically, in the reduced model (Equations 20–22), V and y are fast variables, and z is a slow variable.

4.2. Current-Voltage Relations in the Reduced Model

In this section, by treating the slow variable z as a parameter, we analyze the stability of fixed points with the current-voltage (I-V) relation in the reduced model. In the parameter space of the biophysical condition, we obtain y(t) = V(t) by solving dy/dt = 0 in Equation (21). Substituting it into Equation (20), we get the I-V relation I(V, z) (the full equation, please refer to Supplementary Material: Current-voltage relations). In a voltage-clamp experiment, the I(V, z) measures under the given value of z the membrane current needs to inject to clamp the potential to the value of V. If I(V, z) = 0, the neuron is in the rest or equilibrium state at the value of V. The I-V curve of the reduced model is plotted in Figure 3 for various values of z and I_syn. The equilibrium points are only plotted for z = −65 mV because other curves with different values of z have similar fixed points. The stability of each equilibrium point is evaluated by the linear stability analysis (refer to Supplementary Material: Linear stability analysis).

FIGURE 3

Figure 3. The current-voltage relation in the reduced model. In each panel, we plot multiple I-V curves with different values of z. I-V relations under the different settings of I_syn are plotted in (A–C) separately. The colored points indicate the fixed points when z is set to –65 mV. Other fixed points under the different settings of z are similar with these of z = −65 mV, so we omit them in the figure. (A) I-V curves when $I_{syn}=0.\mu \text{A}/\text{c}{\text{m}}^2$. For the high value of z, the model has three equilibrium points. (B) I-V curves when $I_{syn}=-0.05\mu \text{A}/\text{c}{\text{m}}^2$. The injection of inhibitory current results in the lower value of resting potential V. (C) I-V curves when $I_{syn}=0.10\mu \text{A}/\text{c}{\text{m}}^2$. Depolarizing currents shift the I-V curves upward, making the model only have one isolated unstable point associated with a stable limit cycle (corresponds to the action potential). The neuronal parameters used in this study are $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

Figure 3A presents the I-V curves when no external input is applied. It reveals that for a high value of z, there are three equilibrium points: a stable focus (S) at the low membrane potential which corresponds to the resting potential, and two unstable fixed points (U1 and U2) with high potential values. Usually, the unstable focus (U2) with the highest potential value (near –20 mV) is surrounded by a stable limit cycle attractor. This means the system will exhibit bistability for a high value of z. Gradually decreasing z, the stable equilibrium point S and the unstable point U1 get close to each other, coalesce at a saddle-node bifurcation point, and then disappear. This makes the system display rapid action potentials because only the unstable point U2 and its associated limit cycle are left when the value of z is low.

Injection of the external current changes the position of I-V curves relative to the zero line I = 0. Figures 3B,C show the I-V curves of the reduced model when the inhibitory and excitatory currents are applied, respectively. Biophysically, the inhibitory current usually hyper-polarizes the membrane potential to a lower value. This can be interpreted by the I-V curves shown in Figure 3B. The inhibitory current moves the I-V curves downward, making the resting point S left-shifted. Moreover, the phenomenon that the injection of excitatory current depolarizes the neuron to produce action potentials can also be qualitatively understood by the I-V curve shown in Figure 3C. Because depolarizing currents shift the curves upward, eliminating the bistable property of the model and leaving an isolated unstable point surrounded by a stable limit cycle (corresponds to the action potential).

4.3. Rebound Bursting via the “Fold/Homo-Clinic” Bifurcation

The I-V curves above let us inspect the neuronal excitability under the different values of z and I_syn. However, the mechanism of why the model produces rebound bursting and tonic spiking is still not clear. By continually treating the slow variable z as a bifurcation parameter, in the next two sections, we analyze the phase portraits of the fast subsystem, and then interpret the model as the evolution of the phase portraits under the control of the slow variable z.

The behavior sequence in a rebound bursting (Figure 4A) when the TRN neuron receives a hyperpolarizing current is illustrated in Figure 4B. At the resting potential without external input (①), the model exhibits a bistable behavior: a stable node (corresponding to the resting potential), a saddle-node, and an unstable node with a limit cycle. Once a hyperpolarizing current I = −0.06 μA/cm² is applied (②), the value of the stable node starts to decline, resulting in lower values of the membrane potential V and the slow variable z. Further, the withdrawal of such inhibitory current (③) removes the stable equilibrium state, with only an unstable focus left (associated with a limit cycle). This leads the model to produce rapid action potential and gradually increase the slow variable value z. The increasing z once again makes the model bistable (④). However, the model is still in the stable limit cycle state due to the hysteresis. The repetitive spiking stops when z is too large. This happens when z passes the value of –63.25 mV (⑤), the limit cycle becomes a homo-clinic orbit to a saddle, and then disappears. When z increases to the highest value (⑥), the bursting trajectory jumps down to the stable equilibrium corresponding to the resting state (stable node). Therefore, according to the bursting classification schema of Izhikevich (2007), the firing pattern exhibited in the TRN neuron model should be named the “fold/homo-clinic” bursting, because the resting state disappears through the fold bifurcation (② → ③), then the spiking limit cycle state disappears through the saddle homo-clinic orbit bifurcation (⑤).

FIGURE 4

Figure 4. Phase plane analysis during the TRN rebound bursting. (A) The bursting behavior of the three-variable TRN model after receiving an inhibitory current. The external input with the size of −0.03 mS/cm² is given during 50–250 ms. (B) Transitions of phase portraits as the slow variable changes over time. ①: Without the external input, the model is bistable. Under the initial condition, the model is in the resting state (corresponds to the stable node). ②: Applying an inhibitory current makes us get smaller values of resting potential V and slow variable z. ③: Once upon the current is removed, only the unstable focus and its associated limit cycle (action potentials) are left. ④: Repetitive action potentials increase the slow variable z and make the model bistable again. ⑤: Once z is so large that the limit cycle becomes a homo-clinic orbit to a saddle and the repetitive spiking cannot be sustained. ⑥: When z reaches its highest value, the model directly jumps to its global stable node. The neuronal parameters used in this study are $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

To summarize the phase portraits shown in Figure 4B, a bifurcation analysis for the rebound bursting concerning the slow variable z is presented in Figure 5A. For z < −66.54 mV, the model only exhibits stable action potentials. The amplitude of action potential decreases with the smaller value of z. However, for −66.54 mV ≤ z ≤ −63.33 mV, the model exhibits bi-stability due to the emergence of a stable node from the fold bifurcation. The saddle-node denotes the separation of attraction regions between two stable attractors. For z > −63.33 mV, only the stable nodes are exhibited due to the disappearance of the action potential through the saddle homo-clinic orbit bifurcation.

FIGURE 5

Figure 5. The bifurcation diagram of the TRN neuron. The colored dots indicate the fixed points evaluated in the V-y two-dimensional sub-system (refer to Figure 4B). They constitute the slow manifold. The “x” dots, intersections of the slow manifold and the z nullcline, denote the global fixed points evaluated in the V-y-z system. (A) The bifurcation structure when $I_{syn}=0.0\mu \text{A}/\text{c}{\text{m}}^2$. When z is small, the model produces a barrage of action potentials, i.e., bursting. Meanwhile, z increases until it reache z = –63.33 mV, where a saddle homo-clinic orbit bifurcation occurs. Then z decreases due to dz/dt < 0, and the trajectory terminates due to the attraction of the global stable focus. (B) The bifurcation structure when $I_{syn}=0.1\mu \text{A}/\text{c}{\text{m}}^2$. The depolarizing current changes the bifurcation structure, and results in two distinct manifolds. The first is a region where an unstable focus is surrounded by a stable limit cycle. The second is a regime in which multiple fixed points coexist. During the first regime, the decrease of z (dz/dt < 0) exactly balances the increase of z (dz/dt > 0), leading the model to have the ability to produce tonic spiking. (C) The bifurcation structure when $I_{syn}=-0.03\mu \text{A}/\text{c}{\text{m}}^2$. Under the weak input, the model to have the ability to produce periodic bursts. (D) The bifurcation structure when $I_{syn}=-0.06\mu \text{A}/\text{c}{\text{m}}^2$. When the hyper-polarization current is strong, the model will jump to the stable focus state after a homo-clinic bifurcation. The parameters used in this figure are $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

In Figure 5A, we also plot the z nullcline V(t) = z(t), which is obtained by solving dz/dt = 0 in Equation (22). The intersections of the z nullcline and the bifurcation structure, denoted by the “ × ” symbols, are evaluated by the linear stability analysis (refer to Supplementary Material: Linear stability analysis). When initialized in the attraction domain of limit cycles, the model produces a barrage of action potentials (refer to the “trajectory” line in Figure 5A). During the spike train of action potentials, z increases until it reaches a point (z = −63.33 mV) that corresponds to the saddle homo-clinic orbit bifurcation. Then, the trajectory is attracted toward the stable nodes (red dots in Figure 5A), and z decreases due to dz/dt < 0. z decreases until it goes across the z nullcline, and V starts to increase when a fold bifurcation occurs at z = −66.54 mV. Finally, the trajectory will terminate at the lowest intersection, which is a stable focus point.

4.4. Tonic Spiking Due to the Fold Cycle Bifurcation

Another key feature of the TRN neuron is the tonic spiking (Figure 6A) during a long-lasting excitation. Figure 6B shows the corresponding phase portraits. Same as Figure 4B, in the resting state (without external input) the model has a bistable behavior, in which a stable node and a stable limit cycle coexist (①). Then, the injection of an excitatory current (②) destroys this bi-stability property, making the model only have a stable action potential attractor. Since the slow variable z, which is initially small, gradually increases, the model exhibits the decelerated burst discharges at the beginning of the current injection (②). Once z is big enough, the slow variable z will show periodic orbit with the same frequency as the fast variable V (③ and ④). This stable periodic solution, or the tonic spiking, will last until the external excitatory current is removed (⑤). Because the removal of the excitatory current will change the distribution of fixed points in the model, leaving a global stable node corresponding to the resting state (⑤).

FIGURE 6

Figure 6. Phase plane analysis of the TRN tonic spiking. (A) The model behavior when receiving a constant excitatory input. The external current with the size of 0.1 mS/cm² is injected during 50–350 ms. The model first shows transient decelerated burst and then displays tonic spikes. (B) The evolution of phase portraits as the slow variable changes in time. ①: Without the external input, the model exhibits bistable behavior. Under the initial condition, the model is in the resting state and has a small value of slow variable z. ②: The excitatory input switches the bistable state to an unstable focus associated with a stable limit cycle. Due to the initial slow variable, z is small, the model first exhibits a barrage of fast action potentials. ③ and ④: Once variable z reaches to big value, the potential V and the slow variable z co-evolve with each other, and the model exhibits stable action potentials. ⑤ When the input is withdrawn, a global stable state emerges and the model trajectory jumps to the resting state. The parameters used in this study are $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

Figure 5B summarizes the phase portraits during a TRN tonic spiking. Different from the bifurcation diagram of a rebound bursting (Figure 5A), when an excitatory current with the size of 0.1 μA/cm² is injected, the model displays two distinct manifolds. The first is a region in which the model has an unstable focus surrounded by a stable action potential limit cycle when z is small. This manifold region terminates after a fold limit cycle bifurcation occurs at z = −56.26 mV. Then the second region appears in which multiple fixed points coexist. In this manifold, the stable limit cycle disappears, only leaving a stable node corresponding to the resting potential.

4.5. One-Parameter Bifurcation Analysis

In the above, we have investigated the bifurcation mechanism of the rebound bursting and the tonic spiking by a zero current and an excitatory current, respectively. In the next, we continue to inspect the neuronal excitability by varying the size of external currents.

Figure 5C presents the bifurcation diagram of the TRN neuron model when given an inhibitory current with the size of −0.03 μA/cm². As before, z will increase until z = −63.77 mV, whereupon the action potentials disappear through a saddle homo-clinic orbit bifurcation. z will continue to increase until the trajectory passes through the z nullcline. Then the system moves toward the stable node solution, and z decreases due to dz/dt < 0 until the stable node disappears when the saddle-node bifurcation occurs at z = −69.775 mV. Different from the bifurcation diagram shown in Figure 5A, the lowest intersection of the z nullcline and the slow manifold is an unstable focus, making the system moves away to the attraction domain of limit cycles again. The critical feature in Figure 5C is the bi-stability between two stable solutions (the limit cycle and the stable node) for z between two bifurcations (the saddle-node bifurcation and the saddle homo-clinic orbit bifurcation). It is this bi-stability that is important for the periodic bursting behavior shown in Figure 5C. Moreover, the length of the bistable interval determines the number of spikes in a burst.

However, with the increasing value of the inhibitory current, the unstable intersection becomes stable. In Figure 5D, we present the bifurcation diagram when the current of −0.06 μA/cm² is injected. This time z increases until z = −64.245 mV, and the fold bifurcation occurs at z = −74.8 mV. Once the limit cycle disappears through the homo-clinic bifurcation, the trajectory will directly jump to the steady-state where the z nullcline intersects the slow manifold at −74.776 mV.

A one-parameter bifurcation diagram summarizing the change in the model dynamics with the external input I_syn is shown in Figure 7. For ${\text{I}}_{\text{syn}}<-0.052\mu \text{A}/\text{c}{\text{m}}^2$, the system has a steady-state as previously demonstrated in Figure 5D. As I_syn increases, a super-critical Andronov-Hopf bifurcation (SHB) occurs at ${\text{I}}_{\text{syn}}=-0.052\mu \text{A}/\text{c}{\text{m}}^2$ (point “SHB1” in Figure 7) due to a stable equilibrium state changes to a stable periodic orbit along with a two-dimensional unstable manifold. The stable orbit exhibited in the model is the periodic bursting (the red and green dotted lines ⋯ in Figure 7). This corresponds to the situation in Figure 5C, where the bi-stability between two bifurcations leads the system to have a periodic bursting solution. The number of the spike in a burst is determined by the size of the inhibitory current. Specifically, the smaller size of the injected inhibitory input, fewer spikes in a burst (refer to ① and ② in Figures 7, 8). Intuitively understanding, this is because the higher value of inhibitory current can drive the neuron into a more hyperpolarized state. The more the neuron is hyperpolarized, the low-threshold calcium channel I_T is more deinactivated, further resulting in the neuron producing stronger bursting. Further increasing I_syn makes the model stable again via another SHB at ${\text{I}}_{\text{syn}}=-0.003\mu \text{A}/\text{c}{\text{m}}^2$ (point “SHB2” in Figure 7). The post-inhibitory excitation current with the size in this parameter region (specifically, $-0.003\mu \text{A}/\text{c}{\text{m}}^2<{\text{I}}_{\text{syn}}<0.059\mu \text{A}/\text{c}{\text{m}}^2$) will result in a rebound bursting as illustrated in Figures 4, 5A. 0.059 μA/cm² is the minimum excitation current necessary to trigger an action potential of the model, because the stable limit cycle occurs after a fold cycle bifurcation appears at ${\text{I}}_{\text{syn}}=0.059\mu \text{A}/\text{c}{\text{m}}^2$ (point “FCB” in Figure 7). This corresponds to the situation shown in Figures 5B, 6 where the increase of z during the trajectory is above the z nullcline exactly balances the decrease of z during the trajectory is below the z nullcline, making the sustained action potential occurs (the solid lines—in Figure 7). It is worth noting that during this parameter region (specifically $0.059\mu \text{A}/\text{c}{\text{m}}^2<{\text{I}}_{\text{syn}}<0.1316\mu \text{A}/\text{c}{\text{m}}^2$) the system is bistable, as a stable limit cycle and a stable-steady state coexist. This means the final state of the system depends on the initial conditions (refer to ③ in Figures 7, 8). The bistable mode continues until ${\text{I}}_{\text{syn}}=0.1316\mu \text{A}/\text{c}{\text{m}}^2$, where a fold bifurcation occurs (“FC” point in Figure 7). While the model still exhibits sustained action potentials (refer to example of ④ shown in Figures 7, 8).

FIGURE 7

Figure 7. One parameter bifurcation diagram of the TRN neuron model. “SHB” denotes the super-critical Andronov-Hopf bifurcation, “FCB” presents the fold cycle bifurcation, and “FB” indicates the fold bifurcation. The dotted line is the maximum and minimum values of V during the periodic burst, and the solid line (—) indicates the maximum and minimum values of V during the limit cycle. Every colored dot is evaluated by the linear stability analysis of the V-y-z system. For the detalied explanation please refer to the main text. The parameters used are $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

FIGURE 8

Figure 8. The simulation results under the different settings of I_syn and the initial state. ① and ② illustrate that the number of spikes in a burst is determined by the size of the inhibitory current. ①: when $I_{syn}=-0.03\mu \text{A}/\text{c}{\text{m}}^2$, the model produces the robust periodic bursting in which each burst has three spikes. ②: decreasing the inhibitory current to $I_{syn}=-0.01\mu \text{A}/\text{c}{\text{m}}^2$, the number of spikes in a burst decreases to one. ③ demonstrates the model has the bistable property when depolarizing current I_syn is small. One is the stable resting state, the other is the stable limit cycle. Specifically, when $I_{syn}=0.07\mu \text{A}/\text{c}{\text{m}}^2$, models with different initial states are evolved to different attractors. ④ illustrates that once the external input is big enough, i.e., $I_{syn}=0.14\mu \text{A}/\text{c}{\text{m}}^2$, the model only has a stable limit cycle attractor, which corresponds to the tonic spiking mode as seen in the TRN neuron. The parameters used in this figure are $g_{KL}=0.0065\text{ mS}/\text{c}{\text{m}}^2$, $g_T=2.25\text{ mS}/\text{c}{\text{m}}^2$, and k = 0.01.

5. Conclusion and Discussion

The thalamic reticular nucleus is a strategical locus due to its pacemaking role in spindle generation during sleep (Steriade et al., 1985, 1987) and its gatekeeper function in selective attention during wakefulness (Halassa et al., 2014; Wimmer et al., 2015; Nakajima et al., 2019). Understanding its intrinsic dynamics is an important step toward understanding its functions. Previously, TRN neurons were modeled using complex realistic conductance-based models whose dynamics are hard to analyze and visualize. In this study, we developed a reduced three-variable neuron model to capture the key dynamical features of TRN neurons. We demonstrated that the reduced model can effectively reproduce the characteristic firing patterns of TRN neurons (Figure 2), but its dynamics are much easier to be analyzed. The bifurcation diagrams illustrate the underlying mechanisms of the two characteristic TRN firing patterns, that is, the rebound bursting is “fold/homo-clinic” bifurcation (Figures 4, 5) and the tonic spiking is the fold cycle bifurcation (Figures 5, 6). Furthermore, one-parameter bifurcation analysis demonstrates that the model displays varying degrees of bursting and tonic discharges when the size of the external current changes (Figures 7, 8). We expect that this study will facilitate our future work to study the complicated dynamics of the TRN network.

It has been proved feasible that neuron models with high-dimensional voltage-gating variables by using the Hodgkin-Huxley schema can be described by a reduced system with fewer essential variables (Chay, 1985; Av-Ron et al., 1993; Golomb et al., 1993; Chay et al., 1995; Maeda et al., 1998; Izhikevich, 2007). The pioneering work is the FitzHugh-Nagumo model (FitzHugh, 1961; Nagumo et al., 1962; Izhikevich and FitzHugh, 2006). It provides a simplified model for the Hodgkin-Huxley model (Hodgkin and Huxley, 1952) and allows us to inspect the mechanism of neuronal excitability and spike generation from the geometrical viewpoints. Later, by the massive numerical simulation, Krinskii and Kokoz (1973) provided an empirical conclusion that gating variables n and h in the Hodgkin-Huxley model have a linear relationship: n(t) + h(t) ≈ 0.84. Thus, h and n variables can be approximated by a single variable. Inspired by this reduction idea, Abbott and Kepler (1990) and Kepler et al. (1992) suggested a more general method to reduce the Hodgkin-Huxley-type neuron models. The core concept behind this method is that: if gating variables behave in a similar time scale, a single variable can be obtained by the weighted average of them. Specifically, the weights should be the relative contributions to the overall membrane potential change. In this study, we applied this method to reduce a full TRN neuron model (Bazhenov et al., 1998, 2002) into a three-variable model. Discharge patterns are essential factors in neural information processing. Our reduction model can reproduce the firing patterns seen in the original model. This demonstrated our reduction is valid (Golomb et al., 1993; Maeda et al., 1998).

Since the seminal paper of Rinzel and Ermentrout (1998) introduced a geometrical method for phase plane analysis on neuronal models, phase plane analysis has become standard on neural modeling. Afterward, under the guidance of such dynamical system theory, simplified but powerful neuron models were proposed, such as the Izhikevich neuron model (Izhikevich, 2003, 2007) and adaptive exponential integrate-and-fire model (Brette and Gerstner, 2005; Gerstner and Brette, 2009). These low-dimensional abstract neuron models are capable of producing tonic spiking and bursting patterns. However, the neuronal parameters largely differ between firing patterns, making the continuous switch between different firing patterns (as seen in biological experiments Bal and McCormick, 1993; Llinás and Steriade, 2006; Halassa and Acsády, 2016) difficult. On the contrary, our proposed three-variable model which is reduced directly from the realistic conductance neuron model can effectively capture different firing patterns and the coherent switch between them (Figures 7, 8). The reduced three-variable model retains the fundamental biophysical properties of the original, as each channel dynamics can be easily recovered from the reduced variables [refer to Eq. (S22) in Supplemental Material].

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. All datasets generated for this study are included in the article. The source code for the simulation and analysis is available at https://github.com/chaoming0625/TRNNeuronAnalysis.

Author Contributions

CW and SW designed the project. CW, SL, and SW performed the formal analysis and wrote the manuscript. CW carried out simulations. SW supervised the project. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by Guangdong Province with grant (no. 2018B030338001, SW) and Huawei Technology Co., Ltd. (no. YBN2019105137, SW).

Conflict of Interest

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.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2021.764153/full#supplementary-material

References

Abbott, L., and Kepler, T. B. (1990). “Model neurons: from hodgkin-huxley to hopfield,” in Statistical Mechanics of Neural Networks, ed L. Garrido (Berlin; Heidelberg: Springer), 5–18.

Google Scholar

Ahrens, S., Jaramillo, S., Yu, K., Ghosh, S., Hwang, G.-R., Paik, R., et al. (2015). Erbb4 regulation of a thalamic reticular nucleus circuit for sensory selection. Nat. Neurosci. 18, 104–111. doi: 10.1038/nn.3897

PubMed Abstract | CrossRef Full Text | Google Scholar

Av-Ron, E., Parnas, H., and Segel, L. A. (1993). A basic biophysical model for bursting neurons. Biol. Cybern. 69, 87–95. doi: 10.1007/BF00201411

PubMed Abstract | CrossRef Full Text | Google Scholar

Bal, T., and McCormick, D. A. (1993). Mechanisms of oscillatory activity in guinea-pig nucleus reticularis thalami in vitro: a mammalian pacemaker. J. Physiol. 468, 669–691. doi: 10.1113/jphysiol.1993.sp019794

PubMed Abstract | CrossRef Full Text | Google Scholar

Barthó, P., Slézia, A., Mátyás, F., Faradzs-Zade, L., Ulbert, I., Harris, K. D., et al. (2014). Ongoing network state controls the length of sleep spindles via inhibitory activity. Neuron 82, 1367–1379. doi: 10.1016/j.neuron.2014.04.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazhenov, M., Timofeev, I., Steriade, M., and Sejnowski, T. J. (1998). Cellular and network models for intrathalamic augmenting responses during 10-hz stimulation. J. Neurophysiol. 79, 2730–2748. doi: 10.1152/jn.1998.79.5.2730

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazhenov, M., Timofeev, I., Steriade, M., and Sejnowski, T. J. (2002). Model of thalamocortical slow-wave sleep oscillations and transitions to activated states. J. Neurosci. 22, 8691–8704. doi: 10.1523/JNEUROSCI.22-19-08691.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Beierlein, M. (2014). Synaptic mechanisms underlying cholinergic control of thalamic reticular nucleus neurons. J. Physiol. 592, 4137–4145. doi: 10.1113/jphysiol.2014.277376

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunton, J., and Charpak, S. (1997). Heterogeneity of cell firing properties and opioid sensitivity in the thalamic reticular nucleus. Neuroscience 78, 303–307.

PubMed Abstract | Google Scholar

Chay, T. R. (1985). Chaos in a three-variable model of an excitable cell. Physica D 16, 233–242. doi: 10.1016/0167-2789(85)90060-0

CrossRef Full Text | Google Scholar

Chay, T. R., Fan, Y. S., and Lee, Y. S. (1995). Bursting, spiking, chaos, fractals, and universality in biological rhythms. Int. J. Bifurcat. Chaos 5, 595–635. doi: 10.1142/S0218127495000491

CrossRef Full Text | Google Scholar

Clemente-Perez, A., Makinson, S. R., Higashikubo, B., Brovarney, S., Cho, F. S., Urry, A., et al. (2017). Distinct thalamic reticular cell types differentially modulate normal and pathological cortical rhythms. Cell. Rep. 19, 2130–2142. doi: 10.1016/j.celrep.2017.05.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Contreras, D., Curro Dossi, R., and Steriade, M. (1992). Bursting and tonic discharges in two classes of reticular thalamic neurons. J. Neurophysiol. 68, 973–977. doi: 10.1152/jn.1992.68.3.973

PubMed Abstract | CrossRef Full Text | Google Scholar

Coulon, P., and Landisman, C. E. (2017). The potential role of gap junctional plasticity in the regulation of state. Neuron 93, 1275–1295. doi: 10.1016/j.neuron.2017.02.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Crabtree, J. W. (2018). Functional diversity of thalamic reticular subnetworks. Front. Syst. Neurosci. 12:41. doi: 10.3389/fnsys.2018.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

Crick, F. (1984). Function of the thalamic reticular complex: the searchlight hypothesis. Proc. Natl. Acad. Sci. U.S.A. 81, 4586–4590. doi: 10.1073/pnas.81.14.4586

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Bal, T., McCormick, D. A., and Sejnowski, T. J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070. doi: 10.1152/jn.1996.76.3.2049

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, D., Wang, Q., Su, J., and Xi, H. (2017). Stimulus-induced transitions between spike-wave discharges and spindles with the modulation of thalamic reticular nucleus. J. Comput. Neurosci. 43, 203–225. doi: 10.1007/s10827-017-0658-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez, L. M., Vantomme, G., Osorio-Forero, A., Cardis, R., Béard, E., and Lüthi, A. (2018). Thalamic reticular control of local sleep in mouse sensory cortex. Elife 7:e39111. doi: 10.7554/eLife.39111

PubMed Abstract | CrossRef Full Text | Google Scholar

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi: 10.1016/S0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstner, W., and Brette, R. (2009). Adaptive exponential integrate-and-fire model. Scholarpedia 4, 8427. doi: 10.4249/scholarpedia.8427

CrossRef Full Text | Google Scholar

Golomb, D., Guckenheimer, J., and Gueron, S. (1993). Reduction of a channel-based model for a stomatogastric ganglion lp neuron. Biol. Cybern. 69, 129–137. doi: 10.1007/BF00226196

CrossRef Full Text | Google Scholar

Govindaiah, G., Wang, T., Gillette, M. U., Crandall, S. R., and Cox, C. L. (2010). Regulation of inhibitory synapses by presynaptic d4 dopamine receptors in thalamus. J. Neurophysiol. 104, 2757–2765. doi: 10.1152/jn.00361.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Halassa, M. M., and Acsády, L. (2016). Thalamic inhibition: diverse sources, diverse scales. Trends Neurosci. 39, 680–693. doi: 10.1016/j.tins.2016.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Halassa, M. M., Chen, Z., Wimmer, R. D., Brunetti, P. M., Zhao, S., Zikopoulos, B., et al. (2014). State-dependent architecture of thalamic reticular subnetworks. Cell 158, 808–821. doi: 10.1016/j.cell.2014.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Halassa, M. M., Siegle, J. H., Ritt, J. T., Ting, J. T., Feng, G., and Moore, C. I. (2011). Selective optical drive of thalamic reticular nucleus generates thalamic bursts and cortical spindles. Nat. Neurosci. 14, 1118–1120. doi: 10.1038/nn.2880

PubMed Abstract | CrossRef Full Text | Google Scholar

Herd, M. B., Brown, A. R., Lambert, J. J., and Belelli, D. (2013). Extrasynaptic gabaa receptors couple presynaptic activity to postsynaptic inhibition in the somatosensory thalamus. J. Neurosci. 33, 14850–14868. doi: 10.1523/JNEUROSCI.1174-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrera, C. G., Cadavieco, M. C., Jego, S., Ponomarenko, A., Korotkova, T., and Adamantidis, A. (2016). Hypothalamic feedforward inhibition of thalamocortical network controls arousal and consciousness. Nat. Neurosci. 19, 290–298. doi: 10.1038/nn.4209

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Huguenard, J., and Prince, D. (1992). A novel t-type current underlies prolonged ca (2+)-dependent burst firing in gabaergic neurons of rat thalamic reticular nucleus. J. Neurosci. 12, 3804–3817. doi: 10.1523/JNEUROSCI.12-10-03804.1992

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience. Cambridge, MA: MIT Press.

Google Scholar

Izhikevich, E. M., and FitzHugh, R. (2006). Fitzhugh-nagumo model. Scholarpedia 1, 1349. doi: 10.4249/scholarpedia.1349

CrossRef Full Text | Google Scholar

Jourdain, A., Semba, K., and Fibiger, H. C. (1989). Basal forebrain and mesopontine tegmental projections to the reticular thalamic nucleus: an axonal collateralization and immunohistochemical study in the rat. Brain Res. 505, 55–65. doi: 10.1016/0006-8993(89)90115-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kepler, T. B., Abbott, L., and Marder, E. (1992). Reduction of conductance-based neuron models. Biol. Cybern. 66, 381–387. doi: 10.1007/BF00197717

PubMed Abstract | CrossRef Full Text | Google Scholar

Krinskii, V., and Kokoz, Y. M. (1973). Analysis of equations of excitable membranes-i. reduction of the hodgkin-huxley equations to a second order system. Biophysics 18, 533–539.

Google Scholar

Lee, S.-H., Govindaiah, G., and Cox, C. L. (2007). Heterogeneity of firing properties among rat thalamic reticular nucleus neurons. J. Physiol. 582, 195–208. doi: 10.1113/jphysiol.2007.134254

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, L. D., Voigts, J., Flores, F. J., Schmitt, L. I., Wilson, M. A., Halassa, M. M., et al. (2015). Thalamic reticular nucleus induces fast and local modulation of arousal state. Elife 4:e08760. doi: 10.7554/eLife.08760

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Lopez-Huerta, V. G., Adiconis, X., Levandowski, K., Choi, S., Simmons, S. K., et al. (2020). Distinct subnetworks of the thalamic reticular nucleus. Nature 583, 819–824. doi: 10.1038/s41586-020-2504-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Llinás, R. R., and Steriade, M. (2006). Bursting of thalamic neurons and states of vigilance. J. Neurophysiol. 95, 3297–3308. doi: 10.1152/jn.00166.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeda, Y., Pakdaman, K., Nomura, T., Doi, S., and Sato, S. (1998). Reduction of a model for an onchidium pacemaker neuron. Biol. Cybern. 78, 265–276. doi: 10.1007/s004220050432

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Garcia, R. I., Voelcker, B., Zaltsman, J. B., Patrick, S. L., Stevens, T. R., Connors, B. W., et al. (2020). Two dynamically distinct circuits drive inhibition in the sensory thalamus. Nature 583, 813–818. doi: 10.1038/s41586-020-2512-5

PubMed Abstract | CrossRef Full Text | Google Scholar

McAlonan, K., Cavanaugh, J., and Wurtz, R. H. (2006). Attentional modulation of thalamic reticular neurons. J. Neurosci. 26, 4444–4450. doi: 10.1523/JNEUROSCI.5602-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

McAlonan, K., Cavanaugh, J., and Wurtz, R. H. (2008). Guarding the gateway to cortex with attention in visual thalamus. Nature 456, 391–394. doi: 10.1038/nature07382

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A. (1992). Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. Progr. Neurobiol. 39, 337–388. doi: 10.1016/0301-0082(92)90012-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. the IRE 50, 2061–2070. doi: 10.1109/JRPROC.1962.288235

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakajima, M., Schmitt, L. I., and Halassa, M. M. (2019). Prefrontal cortex regulates sensory filtering through a basal ganglia-to-thalamus pathway. Neuron 103, 445–458. doi: 10.1016/j.neuron.2019.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Pellegrini, C., Lecci, S., Lüthi, A., and Astori, S. (2016). Suppression of sleep spindle rhythmogenesis in mice with deletion of cav3. 2 and cav3. 3 t-type ca2+ channels. Sleep 39, 875–885. doi: 10.5665/sleep.5646

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinault, D. (2004). The thalamic reticular nucleus: structure, function and concept. Brain Res. Rev. 46, 1–31. doi: 10.1016/j.brainresrev.2004.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Rinzel, J. (1985). “Bursting oscillations in an excitable membrane model,” in Ordinary and Partial Differential Equations, eds B. D. Sleeman, R. J. Jarvis (Berlin; Heidelberg: Springer), 304–316.

Google Scholar

Rinzel, J., and Ermentrout, G. B. (1998). Analysis of neural excitability and oscillations. Methods Neuronal Model. 2, 251–292.

Google Scholar

Rinzel, J., and Lee, Y. S. (1987). Dissection of a model for neuronal parabolic bursting. J. Math. Biol. 25, 653–675. doi: 10.1007/BF00275501

PubMed Abstract | CrossRef Full Text | Google Scholar

Rovó, Z., Mátyás, F., Barthó, P., Slézia, A., Lecci, S., Pellegrini, C., et al. (2014). Phasic, nonsynaptic gaba-a receptor-mediated inhibition entrains thalamocortical oscillations. J. Neurosci. 34, 7137–7147. doi: 10.1523/JNEUROSCI.4386-13.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M., and Guillery, R. W. (2001). Exploring the Thalamus. San Diego, CA: Elsevier.

Google Scholar

Steriade, M., Deschenes, M., Domich, L., and Mulle, C. (1985). Abolition of spindle oscillations in thalamic neurons disconnected from nucleus reticularis thalami. J. Neurophysiol. 54, 1473–1497. doi: 10.1152/jn.1985.54.6.1473

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Domich, L., and Oakson, G. (1986). Reticularis thalami neurons revisited: activity changes during shifts in states of vigilance. J. Neurosci. 6, 68–81. doi: 10.1523/JNEUROSCI.06-01-00068.1986

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Domich, L., Oakson, G., and Deschenes, M. (1987). The deafferented reticular thalamic nucleus generates spindle rhythmicity. J. Neurophysiol. 57, 260–273. doi: 10.1152/jn.1987.57.1.260

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., McCormick, D. A., and Sejnowski, T. J. (1993). Thalamocortical oscillations in the sleeping and aroused brain. Science 262, 679–685. doi: 10.1126/science.8235588

PubMed Abstract | CrossRef Full Text | Google Scholar

Swanson, L. W., Sporns, O., and Hahn, J. D. (2019). The network organization of rat intrathalamic macroconnections and a comparison with other forebrain divisions. Proc Nat Acad Sci 116, 13661–13669. doi: 10.1073/pnas.1905961116

PubMed Abstract | CrossRef Full Text | Google Scholar

Vantomme, G., Osorio-Forero, A., Lüthi, A., and Fernandez, L. M. (2019). Regulation of local sleep by the thalamic reticular nucleus. Front. Neurosci. 13:576. doi: 10.3389/fnins.2019.00576

PubMed Abstract | CrossRef Full Text | Google Scholar

Vijayan, S., and Kopell, N. J. (2012). Thalamic model of awake alpha oscillations and implications for stimulus processing. Proc. Natl. Acad. Sci. U.S.A. 109, 18553–18558. doi: 10.1073/pnas.1215385109

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Jiang, Y., Liu, X., Lin, X., Zou, X., Ji, Z., et al. (2021). A justi-in-time compilation approach for neural dynamics simulation in International Conference on Neural Information Processing (Springer).

Wells, M. F., Wimmer, R. D., Schmitt, L. I., Feng, G., and Halassa, M. M. (2016). Thalamic reticular impairment underlies attention deficit in ptchd1 y/- mice. Nature 532, 58–63. doi: 10.1038/nature17427

PubMed Abstract | CrossRef Full Text | Google Scholar

Wimmer, R. D., Schmitt, L. I., Davidson, T. J., Nakajima, M., Deisseroth, K., and Halassa, M. M. (2015). Thalamic control of sensory selection in divided attention. Nature 526, 705–709. doi: 10.1038/nature15398

PubMed Abstract | CrossRef Full Text | Google Scholar