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.


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(McAlonan et al., , 2008Halassa 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. (1985Steriade et al. ( , 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.
To reproduce these characteristic firing patterns, detailed conductance-based neuron models were proposed for TRN neurons (Destexhe et al., 1996;Bazhenov et al., 1998Bazhenov et al., , 2002Vijayan 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.

A FULL MODEL OF TRN NEURONS
The full mathematical model of TRN neurons we consider was proposed by Bazhenov et al. (1998Bazhenov et al. ( , 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: where V is the membrane potential, C = 1 µF/cm 2 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 −4 cm 2 (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 2+ 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 The spike adjusting threshold V Na th = V K th = −55 mV and the calcium channel threshold V T th = −3 mV.
norepinephrine (McCormick, 1992). The ion current term is written as, where the leakage current is modeled as I L = g L (V −E L ), with the leakage channel conductance g L = 0.06 mS/cm 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, 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 mS/cm 2 , g K = 10 mS/cm 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, where θ ∞ is the voltage-dependent steady state, τ θ is the voltagedependent time constant, and φ θ is a constant. For a gating variable θ ∈ {m, h, n}, θ ∞ and τ θ are expressed by the voltagedependent state transition variables α and β, i.e., The equations and values of the above variables are listed in Table 1.

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.

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 θ = θ −1 ∞ (θ ). Using the equivalent potentials, the original full model Equation (1) can be re-expressed as, where From Equation (5), we also obtain

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.
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, where all the coefficients {ρ i }, for i ∈ {V, m, h, n, p, q}, are nonnegative, and they satisfy ρ V + ρ m = 1 and ρ h + ρ n + ρ p = 1.

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 Substituting Equations (6-8) into Equation (12) and approximating the dynamics in the first-order Taylor expansion, we have, where 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, 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, where 0 < k < 1 is an undetermined constant, which we normally set k to be 0.01.

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, 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, 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, The above Equations (20-22)

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).

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 2 , 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 2 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).
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.

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) 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 saddlenode 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).

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 injection of inhibitory current results in the lower value of resting potential V. (C) I-V curves when I syn = 0.10µA/cm 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 mS/cm 2 , g KL = 0.0065 mS/cm 2 , and k = 0.01.
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 ( 1 ), 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 2 is applied ( 2 ), 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 ( 3 ) 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 ( 4 ). 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 ( 5 ), the limit cycle becomes a homoclinic orbit to a saddle, and then disappears. When z increases to the highest value ( 6 ), 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 ( 2 → 3 ), then the spiking limit cycle state disappears through the saddle homo-clinic orbit bifurcation ( 5 ).
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 bistability 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.
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.

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 ( 1 ). Then, the injection of an excitatory current ( 2 ) destroys The external input with the size of −0.03 mS/cm 2 is given during 50-250 ms. (B) Transitions of phase portraits as the slow variable changes over time. 1 : Without the external input, the model is bistable. Under the initial condition, the model is in the resting state (corresponds to the stable node). 2 : Applying an inhibitory current makes us get smaller values of resting potential V and slow variable z. 3 : Once upon the current is removed, only the unstable focus and its associated limit cycle (action potentials) are left. 4 : Repetitive action potentials increase the slow variable z and make the model bistable again. 5 : Once z is so large that the limit cycle becomes a homo-clinic orbit to a saddle and the repetitive spiking cannot be sustained. 6 : 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 mS/cm 2 , g KL = 0.0065 mS/cm 2 , and k = 0.01. 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 ( 2 ). Once z is big enough, the slow variable z will show periodic orbit with the same frequency as the fast 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 µA/cm 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 µA/cm 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 µA/cm 2 . Under the weak input, the model to have the ability to produce periodic bursts. (D) The bifurcation structure when I syn = −0.06 µA/cm 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 mS/cm 2 , g T = 2.25 mS/cm 2 , and k = 0.01. variable V ( 3 and 4 ). This stable periodic solution, or the tonic spiking, will last until the external excitatory current is removed ( 5 ). 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 ( 5 ). 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 2 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. Under the initial condition, the model is in the resting state and has a small value of slow variable z. 2 : 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. 3 and 4 : 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. 5 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 mS/cm 2 , g KL = 0.0065 mS/cm 2 , and k = 0.01.

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 2 . 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 2 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 steadystate 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 I syn < −0.052 µA/cm 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 I syn = −0.052 µA/cm 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 1 and 2 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 I syn = −0.003 µA/cm 2 (point "SHB2" in Figure 7). The post-inhibitory excitation current with the size in this parameter region (specifically, −0.003 µA/cm 2 < I syn < 0.059 µA/cm 2 ) will result in a rebound bursting as illustrated in Figures 4, 5A. 0.059 µA/cm 2 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 I syn = 0.059 µA/cm 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 µA/cm 2 < I syn < 0.1316 µA/cm 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 3 in Figures 7, 8). The bistable mode continues until I syn = 0.1316 µA/cm 2 , where a fold bifurcation occurs ("FC" point in Figure 7). While the model still exhibits sustained action potentials (refer to example of 4 shown in Figures 7, 8).

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(Steriade et al., , 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/homoclinic" 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 highdimensional 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 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 mS/cm 2 , g T = 2.25 mS/cm 2 , and k = 0.01. FIGURE 8 | The simulation results under the different settings of I syn and the initial state. 1 and 2 illustrate that the number of spikes in a burst is determined by the size of the inhibitory current. 1 : when I syn = −0.03 µA/cm 2 , the model produces the robust periodic bursting in which each burst has three spikes. 2 : decreasing the inhibitory current to I syn = −0.01 µA/cm 2 , the number of spikes in a burst decreases to one. 3 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 µA/cm 2 , models with different initial states are evolved to different attractors. 4 illustrates that once the external input is big enough, i.e., I syn = 0.14 µA/cm 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 mS/cm 2 , g T = 2.25 mS/cm 2 , and k = 0.01. 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(Bazhenov et al., , 2002 into a threevariable 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(Izhikevich, , 2007 and adaptive exponential integrate-and-fire model (Brette and Gerstner, 2005;Gerstner and Brette, 2009). These lowdimensional 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.