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

where θ∞′ denotes the derivative of θ_∞ with respect to v_θ.

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 non-negative, and they satisfy ρ_V + ρ_m = 1 and ρ_h + ρ_n + ρ_p = 1.

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 δ_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,

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.

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) form the reduced model of TRN neurons. The details of the model are given in Supplementary Material: The details of 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).

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.

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

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

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.

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.

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

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 Isyn<-0.052 μA/cm2, 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 Isyn=-0.052 μA/cm2 (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 Isyn=-0.003 μA/cm2 (point “SHB2” in Figure 7). The post-inhibitory excitation current with the size in this parameter region (specifically, -0.003 μA/cm2<Isyn<0.059 μA/cm2) 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 Isyn=0.059 μA/cm2 (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/cm2<Isyn<0.1316 μA/cm2) 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 Isyn=0.1316 μA/cm2, 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).