Mechanisms of Dominant Electrophysiological Features of Four Subtypes of Layer 1 Interneurons

Identifying dominant electrophysiological features of L1 interneuron subtypes

Previously, we separated the L1 IN population into canopy cells, NGFCs, α7 cells, and VIP cells based on a combination of genetic markers, morphology, connectivity, and electrophysiological features, such as input resistance, sag, and latency to first spike (Schuman et al., 2019).

Here, we systematically studied electrophysiological features of firing patterns based on the analysis of all the sweeps in our dataset. We first identified four features of interest in this L1 dataset (Fig. 1A). Following the Petilla terminology (Ascoli et al., 2008), these features are IR firing behavior, acceleration of the firing rate (Acc), OB, and Adap of the firing rate. To show how salient these features exist in our dataset, we plot each cell in an e-feature space (Fig. 1B). To compare across cells, we used the sweep where the firing rate is ∼10 Hz during the 100 ms to 1 s of the current step. The IR is represented by the CV ISI (size), the Acc and Adap features are represented by the normalized change of the firing rate (x axis), and the OB is characterized by the onset firing rate (y axis). Although this is a useful tool, it is limited in two ways. First, when the cell fires irregularly, the change in the firing rate may reflect the irregularity but not the trend in the firing rate. Second, the information that resides outside the 10 Hz sweeps is lost in this analysis. Still, different IN subtypes (colors) are separated in this space. Canopy cells show IR in many cells but no other features; NGFCs only show Acc; α7 cells are OB without exceptions; and VIP cells are heterogeneous but mostly show IR and Adap. The traces in Figure 1A are marked by arrows in Figure 1B.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Dominant electrophysiological features of L1 interneuron subtypes. A, Examples of traces that show four different types of features: A1, IR firing, A2, Acc in the firing rate, A3, OB, and A4, Adap in the firing rate. A2, Inset, 1/ISI and the corresponding fitted line. A3, Inset, Zoomed-in traces during the first 100 ms. Traces with ∼10 APs are selected as examples. B, Scatter plot of the L1 INs from four subtypes in an E-feature space (n = 69). Each dot represents one IN. Only the cells with a sweep that has ∼9 APs from 100 ms to 1 s are analyzed in this panel (average firing rate <ν> = nAP/0.9 s = 10 Hz). x axis shows the normalized change of the instantaneous firing rate (ν = 1/ISI) from the last to the first ISI between 100 ms and 1 s. y axis is the onset firing rate, calculated by the first ISI. The size of the point represents the CV ISI (>200 ms) of that sweep. Arrows indicate the cells from A. Dashed boxes represent the “hero” cells to model except the α7 cell (for details, see main text). C, Distribution of the four features in C1, Canopy cells, C2, NGFCs, C3, α7 cells, and C4, VIP cells. *IR2 and Acc2 features are classified by looser criteria than IR and Acc, respectively (see main text). Cells with IR2 or Acc2 are only marked in the canopy cells and NGFCs for simplicity. Eleven of 27 canopy cells, 4 of 21 NGFCs, 9 of 22 α7 cells, and 8 of 13 VIP cells have an IR2 feature. All VIP cells are classified as having an IR or IR2 feature. Five of 27 canopy cells, 8 of 21 NGFCs, 0 of 22 α7 cells, and 0 of 13 VIP cells are classified as having an Acc2 feature. D, Two sweeps from an example cell with an Acc2 feature. The cell shows acceleration in the firing rate when the injection current is low, while Adap when the injection current is high. E, F, **The cells that do not show any feature are indicated as other1 or other2. The other2 cells are late-spiking, but the other1 cells are not.

Next, we set quantitative criteria based on all the sweeps, but not just one sweep, to better determine if each cell shows these four features (Fig. 1C). Here, we include more cells than in Figure 1B, since having a sweep of 10 Hz is not required. Traditionally, irregularity is used to describe a unique, unpredictable firing pattern (Galarreta et al., 2004). This can be further reflected in the abrupt change of the ISI (Fig. 1A1). To quantify it, we define an IR score of a pair of neighboring ISIs as 1 minus the ratio of the smaller ISI and the larger ISI. If one cell has at least three pairs of ISIs during the last 500 ms of the depolarizing pulse that have an IR score bigger than the threshold of 0.4, we consider the cell as having an IR feature. We find the IR score is more robust than the traditional measurement CV ISI in quantifying the IR firing for the following reasons. On the one hand, we have fine recordings around the rheobase, such that we can collect ISI pairs from different sweeps to distinguish the IR firing from pure noise. On the other hand, many L1 INs show Acc or Adap, which is hard to distinguish from IR by CV ISI.

We only use the last 500 ms to avoid the potential false-positive cases because cells with strong Adap or OB can change their firing rate rapidly at the onset of the current injection period. We also define a looser criterion for having an IR2 feature, that is, if a cell does not display an IR feature but has at least one pair of ISIs during the latter 800 ms. If a cell with an IR2 feature may have an IR feature if the recordings were done with finer current steps, but it may not if the ISI pairs with large IR scores are purely introduced by noise. IR is widely observed in L1 INs, except for NGFCs (Fig. 1B,C). Noticeably, all the VIP cells have an IR or IR2 feature.

We consider a cell as having an Acc feature if it shows acceleration in the firing rate in all the sweeps. A looser criterion, which only requires one cell to have at least two sweeps around rheobase to show the acceleration in the firing rate, is applied to consider a cell as having an Acc2 feature. These cells with an Acc2 feature usually show acceleration when the injection current is low, but Adap when the injection current is high (Fig. 1D). Cells with an Acc feature are uniquely observed in NGFCs. Further, most NGFCs are considered as having an Acc or Acc2 feature (Fig. 1C2), and some canopy cells are found to have an Acc2 feature (5 OF 27).

Then, to distinguish OB and Adap, we define OB and Adap based on different time frames during the current step. We consider a cell to have an OB feature if, first, the onset firing rate is >100 Hz (first ISI is lesser than 10 ms) in all the sweeps with at least five APs; and second, the firing rate drops at least 50% in the first 100 ms in all the sweeps. We consider a cell to have an Adap feature if the firing rate drops at least 20% from 100 ms to 1 s in all sweeps. The OB feature is observed in most of the α7 cells (20 of 22, the two α7 cells without OB features are just below the OB cutoff in Fig. 1B), and some of the VIP cells, while Adap is observed in most of the VIP cells and some of the α7 cells (5 of 22, Fig. 1C3,C4).

The cells with no associated features are marked with other1 or other2. While the other1 canopy cells show no signature (Fig. 1D), the other2 NGFCs are all late-spiking (Fig. 1E). The passive parameters are listed in Table 2 and shown in Figure 2. The detailed characterizations of individual cells are listed in Extended Data Table 1-1.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Electrophysiological features of the L1 interneurons. Bars above suggest the corresponding pair is significantly different (p < 0.001, Mann–Whitney U test). Mean ± SD. Scatters shows the values measured from individual cells. Jitters are added for better visualization. Rheobase and FI slope are measured by fitting the spike-count curve to an ReLU function. Lines are aligned by the rheobase. Bold black lines indicate the ReLU function with the mean FI slope. Gray lines indicate individual spike-count curves of individual cells.

In this study, we want to have parsimonious models that reproduce the spectrum of these behaviors and help understand the mechanisms behind these behaviors. To achieve these, instead of tuning our model based on the averaged features of a subtype, we choose one “hero” cell from each subtype that only shows one type of the electrophysiological feature of interest and reproduce the behavior of that cell. However, readers should be aware that our model does not consider heterogeneity within each cell type. The analysis on individual cells is provided in Extended Data Table 1-1 if the heterogeneity itself is of interest. Canopy, NGFC, and VIP “hero” cells are marked by the dashed boxes in Figure 1B. The closest sweep of α7 “hero” cell has 17 APs after 100 ms, with CV ISI = 0.09, dν /<ν> = 0.18, and onset firing rate = 235.3 Hz.

We do this because the averaged “cell” may not be physically plausible because of the large heterogeneity within an interneuron subtype. We choose to make individual neuron models based on the EIF neuron (Fourcaud-Trocmé et al., 2003), but not Hodgkin-Huxley models, to minimize the number of parameters we need to tune. In the following sections, we first reproduce IR in a canopy cell model; second, capture the Acc and delayed spiking in an NGFC model; and third, generate the OB in an α7 cell model. Finally, we model a VIP cell that shows IR, OB, and Adap simultaneously.

IR firing is reproduced in a canopy cell model with an SIK channel

Since canopy cells with an IR feature do not have any other of the features under consideration, modeling IR on canopy cells does not need to consider the potential interactions between different features (Fig. 1C1). From the selected “hero” canopy cell, we observed substantial subthreshold fluctuations ∼–40 mV from the voltage traces just above the rheobase (Fig. 3A,C). This is more obvious in the phase diagram of the cell (Fig. 3B,D), where we plot the deviation of the voltage, which is proportional to the net current into the cell, over the voltage. A quasi-stable point is observed ∼–40 mV where the net current goes <0, suggesting some outward current is activated during this range.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Modeling IR spiking of a canopy cell with an SIK channel. A, The voltage trace of an example canopy cell at I_inj = 376 pA. B, Phase diagram of the cell. Inset, The corresponding zoom-out diagram. y axis is the slope of the voltage trace, which is proportional to the net current to the cell. Red arrows indicate where we estimate the value of reset, curvature, and firing threshold of the corresponding model. C, D, The trace and the phase diagram at I_inj = 380 pA. E, The sketch of the canopy cell model (left) and the equilibrium value of the gating variables of the SIK channel (right). F, Δb estimation. Left, The shapes of all the APs from all the sweeps of the same cell. We exclude the first two APs in every sweep. APs are aligned at the maximum of the voltage. Dashed lines indicate t = –0.4, 1.6 ms, respectively. Voltage traces of these 2 ms time windows are used to calculate the jump of the slow inactivation variable Δb. Right, The histogram of Δb with an initial value b_init = 0.28. The value b_init = 0.28 is the average value after 3 s of injection in the simulation. The average change is <Δb> = –0.0020. G, The voltage SD measured between 500 ms and 1 s during the current injection from the data (black) and the model (red). H, I, The voltage trace and the phase diagram of the canopy cell model at I_inj = 370 pA. J, The voltage trace at I_inj = 380 pA. K, The comparison of irregularity from the data and the model by IR score (left) and by CV ISI (right). L, The comparison of firing rate from the data, estimated by the number of APs during the 1 s current injection, and from the model, estimated by the steady-state firing rate. Red arrows indicate I_inj = 376, 380 pA. M, The initial dynamics are not considered in this model. Top, The voltage trace that includes the initial transient part when I_inj = 380 pA. Bottom, The transient dynamics of gating variables of the SIK channel. Black and red lines indicate the fast and slow variables a and b, respectively.

To mimic the IR at the latter 500 ms during the current injection, we include an SIK channel in an EIF model (Fig. 3E, left). The major parameters of an EIF model, namely, reset voltage, firing threshold, and curvature at the firing threshold, can be estimated from the phase diagram directly (Fig. 3B, red arrows). The slow-inactivation gating variable from the SIK channel can introduce the clustered spiking, which has been known (Wang, 1993). Here, the model of the SIK channel follows the model of a Kv1 channel in PV⁺ interneurons (Sciamanna and Wilson, 2011), which contains one fast activating gating variable a and one slow inactivating gating variable b. We shift the equilibrium value of the gating variables to a more depolarized value such that the window current of this channel is ∼–40 mV (Fig. 3E, right). As argued by Sciamanna and Wilson (2011, their Fig. 8), this is plausible by taking into account that the Kv1 channels are localized on the axon initial segment but not on the soma.

To tune the model, we further consider the change of the gating variable during the APs. We can estimate the change of the slow variable Δb after each AP by solving the differential equation of b while replacing the voltage term by the shape of APs (Fig. 3F). The initial value of b is chosen based on our simulation. This is plausible because the dynamics of b are much slower (time constant τ_b = 150 ms) than the time window of APs (refractory period of the model 2 ms). Choosing different Δb in the model has little impact on the simulation results. Since the jump of the fast variable Δa does not change the dynamics of the model because of the rapid convergence to the equilibrium point (τ_a = 6 ms, Fig. 3M), we set Δa = 0. Also, we estimate the noise term σ based on the SD of voltage before or around the rheobase current during 500 ms to 1 s of the depolarization current step (Fig. 3G). Importantly, we tune the conductance g_SIK of the SIK channel based on the IR score curve. At last, we adjust all the parameters in the model to best reproduce other cell features, namely, the rheobase, the slope of the f-I curve, and the firing behavior of individual sweeps.

The resulting model can reproduce the IR behavior observed in the cell (Fig. 3H–J) and the changing of IR behavior when the injection current varies (Fig. 3K). The firing rate is also comparable between the data and the model (Fig. 3L).

However, since we focus on the irregularity in this section, the dynamics around the beginning of the current injection are not modeled (Fig. 3M). As a result, the model shows a delayed-firing feature (Fig. 3M, top), because of the slow inactivation of SIK channel (Fig. 3M, bottom).

Mechanism of IR firing in the canopy cell model

In our model, The IR feature is reproduced by including an SIK channel. To understand the mechanism, we did a fast-slow analysis on our model (Fig. 4). Excluding the noise term by setting σ = 0, we observe clustered spiking in the model (Fig. 4A). We further assume the slow gating variable b is a constant and analyze the dynamics of the fast manifold defined by the voltage v and the fast activation variable a. As expected for the EIF model, the system has one stable state (SS) and one unstable steady state (US) when the injection current is low or the outward current is large (Fig. 4B). The SS branch and US branch emerge through a saddle-node bifurcation on an invariance cycle (SNIC) by increasing b over the bifurcation point b_min (Fig. 4C). In the simulation of the full system (Fig. 4A,C), the slow variable b changes across the bifurcation point b_min over time. Thus, the full system oscillates between a resting state (the global SS exists) and a spiking state (the global SS vanishes). This dynamic is known as square-wave bursting (Rinzel, 1987; Wang and Rinzel, 1995; Rinzel and Ermentrout, 1998). Adding noise to the system σ > 0, the oscillation between the resting and spiking states mimics the IR firing patterns.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

A fast-slow analysis of the canopy cell model shows that the irregularity results from the square-wave bursting. A, Sample simulation when I_inj = 375 pA without noise. Black line indicates the voltage trace. Red line indicates the difference between b and b_min. Red dashed line is at b – b_min = 0, suggesting the boundary between the resting state and the firing state of the cell. B, Nullclines when I_inj = 375 pA, b = 0.27. The intersections indicate the fixpoints of the fast manifold. C, The bifurcation diagram when varying b with I_inj = 375 pA. The system undergoes an SNIC. Black solid line indicates the stable branch. Dashed line indicates the unstable branch. Blue line indicates the simulation of the canopy cell model while we remove the noise term σ = 0. Red dashed line indicates the minimum b that the fast manifold has a global SS. Arrows indicate the direction of the dynamics. C1, Dynamics of the entire AP range. C2, Zoomed-in dynamics around the bifurcation point.

Spike frequency acceleration and delayed firing are reproduced in an NGFC model

After capturing the dynamics of IR, we next model the acceleration observed in an NGFC. As shown in Figure 5, the Acc is accompanied by the late onset of spiking (Fig. 5B,C). Similarly, as in tuning a canopy cell model, we further characterize the cell based on the corresponding phase diagram (Fig. 5D) where a similar quasi-stable point is observed at ∼–40 mV, the jump of change of the slow variable Δb after each AP (Fig. 5E), and the noise level (Fig. 5F). Importantly, the firing rate slowly accelerates over the 1 s injection period (Fig. 5G,H).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Modeling the acceleration and late-spiking of an NGFC with a smaller SIK conductance. A-C, Voltage traces of the selected NGFC at I_inj = 198, 188, 190 pA. D, The phase diagram of the cell at I_inj = 190 pA. E, The histogram of Δb with an initial value b_init = 0.4. The average change is <Δb> = –0.0040. F, The voltage SDs of the data and the model. G, The instantaneous firing rate ν = 1/ISI over time. We fit the ISIs (dot) to an exponential decay curve (line) for different injection currents I_inj = 188, 190, 194, 198, 202, 206 pA (color). H, The initial firing rate (black) and the last (red) firing rate at different injection currents. I, The sketch of the NGFC model. J-M, Model results correspond to B-D, A. N, The time of the first spike decreases fast when increasing the injection current in both the data and the model. O, P, The instantaneous firing rate ν increases over time in the model. Organized as in G, H.

Since the delayed firing is associated with the SIK current (Storm, 1988), along with our observations, we tune an EIF model with the same SIK channel that has the identical dynamics as in the canopy cell model to reproduce the ACC feature (Fig. 5I). Compared with the canopy cell model, we reduce the SIK conductance g_SIK to avoid the IR. The relationship between SIK conductance g_SIK and the IR will be discussed later in Figure 11. After a similar tuning procedure, the resulting model is capable of reproducing the Acc feature observed in an NGFC (Fig. 5J–P). In addition, the time of the first AP drops fast when it increases the injection current, and our model is capable of reproducing this behavior (Fig. 5N).

The delayed firing and the acceleration come from the slow inactivation gating variable b of the SIK current. To directly show that, we simulate the model with two sequential voltage steps and current steps with a variable resting time Δtime (Fig. 6A,B). During the first voltage step (Fig. 6A), the model neuron is clamped at –40 mV. The fast-activating gating variable a opens first and increases the SIK current. Then, the slow-inactivating variable b gradually closes and decreases the SIK current. The difference between the onset and steady SIK current is ∼15 pA in our NGFC model. During a rheobase current step (Fig. 6B), the model neuron depolarizes to a subthreshold voltage but is held there by the fast-activating SIK current. It is only released to fire when the SIK current inactivates enough. During a current step that is above the rheobase, the SIK current also inactivates gradually, which leads to an acceleration of the firing rate. Further, the strength of the SIK current is subjected to the initial state of the slow variable b. If the second voltage step or current step is close to the first, the slow variable b does not fully restore to the open state. As a result, the peak of the SIK current and the delay of the first spike is smaller in the second step compared with the first step (Fig. 6C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Recovering dynamics of the NGFC model. A, Simulation of the two voltage steps V_clamp = –40 mV with a resting time Δtime = 200 ms. From top to bottom, the voltage of the model, gating variables a (black) and b (red), SIK current I_SIK. B, Simulation of the two current steps I_inj = 188 pA with a resting time Δtime = 200 ms. From top to bottom, injection current, gating variables a (black) and b (red), the voltage of the model. C, The peak SIK current and delay of the first spike recover with a long resting time Δtime. The recover ratio of the peak SIK current (red) is calculated by comparing the maximum I_SIK during the second voltage step to that during the first voltage step. The delay of the first spike is calculated from the current steps (black).

OB is reproduced in an α7 cell model

We then studied the OB behavior, which is observed in most of the α7 cells. Cells with an OB feature show a characteristic high onset firing rate and a fast dropping of the firing rate during the first 40 ms (Fig. 7A). In the phase diagram of the cell, several lines are above the spiking cycle (Fig. 7B, red arrow), which indicates the existence of a transient current that is less dependent on the voltage. Further, the onset firing rate is ∼270 Hz across sweeps with different injection currents. The firing rate robustly drops below 50 Hz within 40 ms (Fig. 7C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Modeling the OB of an α7 cell with a T-type calcium channel. A, The voltage trace at I_inj = 262 pA. B, Phase diagram of the cell. Red arrow indicates the trajectory from the first AP within the OB period. C, Instantaneous firing rate shows a sudden drop in the first 40 ms of the injection period. Colors represent sweeps with different injection currents, I_inj = 252, 262, 272, 282, 292 pA. D, The sketch of the α7 cell model. E, The voltage SDs of the data and the model. F, The firing rates from the data and the model. G-I, Modeling results organized as in A-C.

As previously suggested (Schuman et al., 2019) and demonstrated (Schuman et al., 2021), OB in α7 cells is dependent on the expression of a T-type Ca²⁺ channel. This was also demonstrated in bursting L2/3 VIP cells (Prönneke et al., 2020). We thus include a T-type Ca²⁺ current in our EIF model for an α7 cell model (Fig. 7D) based on the reported dynamics of T-type Ca²⁺ channels from a thalamocortical relay neuron model (Smith et al., 2000). After tuning the T-type Ca²⁺ channel conductance g_T to the OB in the data, along with noise level σ (Fig. 7E) and other passive parameters, we can reproduce the high onset firing rate and the fast dropping of the firing rate during the first 40 ms in the α7 cell model (Fig. 7F–I).

IR, OB, and adap are simultaneously reproduced in a VIP cell model

After reproducing the IR, Acc, and OB features individually, we turn to the complex dynamics observed in VIP cells. Most VIP cells show more than one considered feature (Fig. 1B4). For example, OB (Fig. 8A,E), irregularity (Fig. 8B,G), and Adap (Fig. 8D,H) are observed in the same selected VIP cell. All these features can be simultaneously reproduced with an EIF model with a T-type Ca²⁺ channel, an SIK channel, and a spike-triggered small conductance K⁺ channel (Fig. 8I).

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

Reproducing the OB, IR firing patterns, and adaption simultaneously in a VIP cell model. A-D, Voltage traces at I_inj = 91, 121, 261 pA, which are examples to show OB, IR, and Adap, respectively. E, The instantaneous firing rates ν from 0 to 100 ms that show the OB. Colors represent sweeps with different injection currents I_inj = 91, 131, 171, 212, 251, 291 pA. F, The histogram of Δb with an initial value b_init = 0.25. The average change is <Δb> = –0.0009. G, The ISI ratio (black) and CV ISI (red) curves that show the irregularity. H, The instantaneous firing rates ν between 100 ms and 1 s that show the Adap. Color coding is the same as in D. I, Sketch of the VIP cell model. J-P, Corresponding modeling results. M, SD of the voltage of the data and the model.

As described in the modeling of an α7 cell (Fig. 7), the OB has the characteristic high onset firing rate and a fast dropping in the firing rate during the first 40 ms. This was reproduced by tuning the conductance g_T of the T-type Ca²⁺ channel, as in the α7 model.

In addition, the IR score is highest in the sweeps around rheobase, and it decays slowly with increasing injection current (Fig. 8G). The IR score remains >0.4, which is the threshold to detect IR, even in a 1 s sweep with 60 APs with injection current I_inj = 300 pA. Above, we reproduced IR spiking by including a strong SIK current in the canopy cell model. However, the IR score drops fast in the canopy cell model when increasing the injection current (Fig. 3K).

Further, strong Adap is also observed in the VIP cell in the last 900 ms (Fig. 8D), where the firing rate drops exponentially to a saturation value (Fig. 8H). To quantify that, we define the AI as the fraction of the firing rate that drops during this 900 ms period. For this VIP cell, AI is ∼0.5 in all the sweeps (Fig. 9E). To reproduce the Adap, we include a spike-triggered small conductance K⁺ channel in our model (Fig. 8I) that generates after-hyperpolarization current reproducing current-independent adaptation, as in Liu and Wang (2001).

• Download figure
• Open in new tab
• Download powerpoint

Figure 9.

More results of the VIP cell model. A, B, The voltage traces of the model in 3 s at I_inj = 121, 141 pA. C, The voltage traces of the model at I_inj = 235, 237 pA. The jump in the curve of the initial firing rate is because of the diminishing of the silent period after the OB. D, The comparison of the firing rates ν of the data and the model. The ν of the model is measured as the number of APs during 1 s current injection. The initial and last (steady-state) firing rates ν of the model are measured from 100 to 1500 ms. E, The AI from the data (black) and the model (red), measured based on the fitting of 1/ISI from 100 ms to 1 s. F, Voltage traces when varying the time constant of fast activation variable τ_a = 40, 80 ms with an injection current I_inj = 141 pA. The firing pattern around 100 ms is more comparable between the data and the model. However, τ_a is not biologically plausible anymore (original τ_a = 4 ms).

After a similar tuning of the VIP cell model, based on the phase diagram (Fig. 8C), the noise level (Fig. 8M), the jump of the Δb (Fig. 8F), we can reproduce simultaneously the OB, which is reflected in the high onset and the fast-dropping of the firing rate during the first 40 ms (Fig. 8J,N), IR, which is reflected by the high IR score and high CV of the ISIs (Fig. 8K,O), and Adap, which shows exponentially decay during the latter 900 ms (Fig. 8L,P).

However, there are limitations to this VIP cell model. First, although the irregularity during the steady state is similar between the data and the model, the onset of IR firing in the model is usually later than that in the data. Especially when the injection current is around the rheobase, the IR firing may be observed after 1 s in the model (Fig. 9A). In contrast, in the data, it may be observed between 500 ms and 1 s (compare Fig. 8B). Second, the cell in the model becomes quiescence after the OB period (∼100 ms, Fig. 9B,C), which results from the strong Adap current triggered by several APs during the OB and the fast activation of the SIK current. But this is not the case in the data (Fig. 8B). In the model, this quiescence period diminishes when the injection current is high enough I_inj = 237 pA (Fig. 9C). Since we are measuring the Adap starting at 100 ms, the diminishing of the quiescence period is reflected in the onset firing rate after 100 ms (Fig. 9D) and is further reflected in the jump on the AI curve (Fig. 9E).

To get an insight into what impacts firing after the OB period, we change the timescale of the fast activation variable of the SIK channel, which delays the onset of outward SIK current (Fig. 9F). As expected, the firing after the OB continues. But to achieve a comparable length with the data, the timescale needs to be ∼80 ms, which is unlikely for an SIK channel. We discuss this discrepancy further in the discussion section.

Mechanism of IR firing in the VIP cell model

To investigate the origin of IR in the VIP cell model, we do a similar fast-slow analysis as we did for the canopy cell model. We first exclude the noise term in the model by setting σ = 0 and observe cluster spiking in the simulation (Fig. 10A). In contrast to the canopy cell model, we have two slow variables: slow inactivation variable b and the dimensionless Ca²⁺ concentration C_ca. To do the fast-slow analysis, we assume both are constants. Similarly, the system has one SS and one US when either outward current prevails the injection current (Fig. 10B). By assuming C_ca is a constant, we can define the b_min(C_ca) as the minimal b where a global SS exists (Fig. 10C, red dashed lines). In the full system, the line b_min(C_ca) represents the line where the SNIC bifurcation happens (Fig. 10D, red dashed lines). Only in the region above the bifurcation line, a global SS exists. In the simulation of the full model (Fig. 10A,D), the system oscillates between a resting state and a spiking state, which is represented as cluster spiking. The dynamics change to IR firing when adding noise to the system by setting σ > 0.

• Download figure
• Open in new tab
• Download powerpoint

Figure 10.

The fast-slow analysis of the VIP cell model shows that the irregularity results from the square-wave bursting. A, The sample simulation when I_inj = 170 pA. Black line indicates the voltage trace. Red line indicates the difference between b and b_min. Red dashed line is at b – b_min = 0, suggesting the boundary between the resting state and the firing state. B, Nullclines of the fast manifold when I_inj = 170 pA, C_ca = 10, b = 0.3. The intersections indicate the fixpoints of the fast manifold. C, The bifurcation diagram when I_inj = 170 pA, C_ca = 10. Black solid line indicates the stable branch. Dashed line indicates the unstable branch. Red dashed line indicates the minimum of b that the fast manifold has a global SS. D, The bifurcation diagram of the whole system and the simulation when I_inj = 170 pA. Red dashed line, b_min(C_ca), indicates the SNIC bifurcation line. Only the region above which has a global SS. Black line indicates the simulation of the VIP cell model without noise (σ = 0).

Mechanism behind the traditional alias “IR spiking cell” for VIP cells

Traditionally, the term IR spiking has been used mainly for VIP cells but rarely for other INs (Tremblay et al., 2016). But, in our study, the irregularity is not uniquely found in VIP cells but distributed in multiple subtypes of L1 interneurons. Especially, the IR in the VIP cells is not because of a higher noise level. Indeed, the noise levels around the rheobase are comparable across all the subtypes (Figs. 3G, 5F, 7E, 8M). Is there anything unique for the IR VIP cells but not other INs with an IR feature?

To answer this question, we systematically vary the parameters of our VIP cell model. We measure the IR score while varying Adap strength γ and the SIK conductance g_SIK at different injection currents (Fig. 11A,B). The IR score is higher when the SIK conductance is bigger. Interestingly, the high IR score is more spread out with a stronger Adap. To quantify this effect, we define IR width as the range of injection current where the IR score is bigger than the threshold 0.6 (Fig. 11C). As a result, we observed that the IR width only starts increasing when g_SIK reaches at ∼100 nS (Fig. 11D, top). However, the IR width remains small if an Adap current is not included γ = 0. The IR width only becomes large when the Adap strength is also increased. The results are similar if we use the IR area under the curve instead of IR width (Fig. 11D, bottom).

• Download figure
• Open in new tab
• Download powerpoint

Figure 11.

The IR firing of VIP cells is more easily observed in experiments because of a large IR width introduced by the interaction between I_SIK and I_adap. A, The ISI ratio (color) over different injection currents and Adap strength γ when g_SIK = 150 nS. B, Same as in A, but g_SIK = 50 nS. C, The sample ISI ratio when g_SIK = 150 nS, γ = 0.1. The range of injection current when CV ISI is above the threshold 0.6 is defined as the IR width. D, The irregularity visibility is high only when the g_SIK and γ are both large. Top, The IR width curves over g_SIK. Different lines indicate results with different Adap strength γ. Dots and arrows indicate the parameters of the VIP cell model (top right), the canopy-like case (bottom right, shown in E), and the NGFC-like case (left, shown in F). Bottom, The IR area-under-curves, defined as the area under the curve of ISI ratio and above the ISI ratio threshold 0.6. E, The model results when set γ = 0, g_SIK = 150 nS look like a canopy cell. Left, The IR score. The drop of the curve is sharp when increasing the injection current, and the IR width is small. Right, Voltage trace at I_inj = 121 pA. F, The model results when setting γ = 0, g_SIK = 90 nS looks like an NGFC. The IR width is further reduced compared with A. In this case, we remove the T-type CA²⁺ current for better visualization.

These explain the behavioral differences observed in the canopy, NGFC, and VIP models. In our VIP cell model, the IR width is at ∼120 pA (Fig. 11D, right top arrow; g_SIK = 150 nS, γ = 0.1). When the Adap current is removed by setting γ = 0, the IR width decreases drastically to ∼20 pA (Fig. 11D, right bottom arrow), where the dynamics mimic the IR in canopy cells (Fig. 11E). If, in addition, we decrease SIK conductance to g_SIK = 90 nS, the IR width decreases to ∼10 pA (Fig. 11D, left arrow), where the dynamics mimic the Acc in NGFC (Fig. 11F). These results suggest that IR is most easily observed only when the g_SIK is big and the Adap is strong, which is the case for VIP cells, but not other INs.

To verify this claim in our dataset, we measured the IR width by multiplying the number of sweeps with an IR score bigger than 0.4 by the corresponding current step size of IR cells in different IN subtypes. The average IR width of VIP cells is the highest, and it is much higher than that of canopy cells (Table 2). Furthermore, the input resistance is the highest of the VIP cells compared with other L1 INs (Table 2), including α7 cells. Although the IR width of VIP cells and α7 cells is comparable, the corresponding voltage change of VIP cells is much higher than that of α7 cells, which explains why the IR feature is most salient in VIP cells.

VIP cell model resonant with theta/alpha input through the activation of T-type Ca²⁺ channel

Because of the intrinsic properties, neurons may respond differently to an input of varied frequencies although the mean is the same. Theoretical and experimental studies suggested that these responses, which are referred to as the DG, are important in understanding the response from a neuron ensemble (Fourcaud-Trocmé et al., 2003; Geisler et al., 2005; Köndgen et al., 2008).

Recently, Borda Bossana et al. (2020) reported that rat L1 INs are homogeneous and narrow-tuned through the DG analysis by using a colored noisy input generated by an OU process, developed by Higgs and Spain (2009; see also Ilin et al., 2013). Briefly, in this analysis, a subthreshold colored noisy current signal, generated from an OU process with a time constant τ_σ = 5 ms, is given as input. Since the average input I_inj is subthreshold, the responding APs are results of noise, but not a tonic, component of the input. The average current I_inj and noise level σ are tuned such that the average firing rate is about <ν> = 5 Hz and the SD of voltage is ∼4 mV. These procedures are required to have similar membrane potential fluctuations that recorded in vivo (Ilin et al., 2013). Next, long recordings are made (at least 1000 s long and collecting 5000 APs), and further analyzed in the Fourier space. The ratio between the response and the input at a specific frequency is defined as the DG at that frequency (DG(f), referred to as the dynamic transfer function in Borda Bossana et al., 2020).

Here, we test the response of our VIP cell model, along with other L1 IN models, to the different rhythm inputs through the DG analysis. We find that, with a designed synaptic input fluctuation (tuned at the black arrow in Fig. 12A), the VIP cell shows two non-zero firing rate branches (Fig. 12A). The sub-rheobase branch is expected where the APs are triggered by colored noise. However, when we decrease the average injection current while keeping the noise σ the same, the second branch around the resting voltage emerged, with an increased potential fluctuation and a peak ∼–60 mV. This second branch disappears if the T-type Ca²⁺ channel is removed from the model (g_T = 0, Fig. 12A, gray lines). Within the sub-rheobase regimen (Fig. 12B, top), the DG shows a similar “narrow bandwidth” as reported (Borda Bossana et al., 2020). The results are robust when the average firing rate varies (3 and 7 Hz, Fig. 12B, top, two thin lines). Interestingly, around the resting regimen (Fig. 12B, bottom), the model resonates to the theta/alpha (4-15 Hz) band input. The observed resonance has characteristics that agree with the dynamic of the T-type Ca²⁺ channel. This resonance happens at the average voltage at ∼–60 mV, where the T-type Ca²⁺ channel activates (Fig. 12A, bottom). Also, the alpha/theta rhythm matches the activating time constant of the T-type Ca²⁺ channel τ⁺_h = 100 ms.

• Download figure
• Open in new tab
• Download powerpoint

Figure 12.

The VIP model resonates with theta/alpha band input only when the average input is around the resting. A, Simulation results of the VIP model when varying the average injection current. In this panel, σ = 50 pA. Gray lines model results while excluding the T-type Ca²⁺ channels. Top, The average firing rate ν shows two non-zero branches. Middle, The voltage's SD decreases when the injection current increases. Bottom, The average voltage of the model. B, The normalized DG when the injection current is at the sub-rheobase region (I_inj = 20 pA, black arrow in A) and around the resting (I_inj = 105 pA, red arrow in A). The DG shows an unexpected high-gain range around the theta/alpha band when the model is around the resting regimen. The high-frequency profile was described as f^–α with a cutoff frequency f₀ and α, shown as a dashed gray line. The two thin black lines indicate the simulation with average firing rates <ν> = 3, 7 Hz, respectively. Solid gray line indicates the 95% confidence level generated by bootstrapping the shuffled ISIs. Top, f₀ = 58.9 Hz, α = 0.99. Bottom, f₀ = 67.6 Hz, α = 0.69. C, Simulation results of VIP model with reduced noise σ = 5 pA. Top, The average firing rate. Bottom, The SD of the voltage. D, The DG around the resting still shows a high gain range around the alpha band (I_inj = 25 pA, red arrow in C) when the input noise is drastically reduced (f₀ = 41.7 Hz, α = 1.52). E-G, The normalized DG for the canopy cell model (I_inj = 282 pA, σ = 90 pA, f₀ = 77.6 Hz, α = 0.62), NGFC model (I_inj = 170 pA, σ = 80 pA, f₀ = 14.8 Hz, α = 0.65), and α7 cell model (I_inj = 200 pA, σ = 65 pA, f₀ = 38.9 Hz, α = 0.51) when the average input current is at sub-rheobase. H, The α7 cell model also has a non-zero firing branch around the resting that shows high gain around the theta/alpha band frequency (I_inj = 25 pA, σ = 65 pA, f₀ = 83.2 Hz, α = 1.10).

To make the voltage fluctuation comparable, we reduced the input noise σ, while the non-zero firing branch around the resting remains (Fig. 12C). Similarly, we observed high gain to the theta/alpha band input (Fig. 12D).

We further tested the DG in other IN models (Fig. 12E–H). All model shows a similar “narrow bandwidth” feature at the sub-rheobase regimen, suggesting the DGs of L1 IN models are homogeneous, as reported (Borda Bossana et al., 2020). Only the α7 model, which includes the T-type Ca²⁺ channel, shows another non-zero firing regimen around the resting, where the model resonates to the theta/alpha band input (Fig. 12H).