## Abstract

Gamma oscillations are believed to play a critical role in in information processing, encoding, and retrieval. Inhibitory interneuronal network gamma (ING) oscillations may arise from a coupled oscillator mechanism in which individual neurons oscillate or from a population oscillator in which individual neurons fire sparsely and stochastically. All ING mechanisms, including the one proposed herein, rely on alternating waves of inhibition and windows of opportunity for spiking. The coupled oscillator model implemented with Wang–Buzsáki model neurons is not sufficiently robust to heterogeneity in excitatory drive, and therefore intrinsic frequency, to account for *in vitro* models of ING. Similarly, in a tightly synchronized regime, the stochastic population oscillator model is often characterized by sparse firing, whereas interneurons both *in vivo* and *in vitro* do not fire sparsely during gamma, but rather on average every other cycle. We substituted so-called resonator neural models, which exhibit class 2 excitability and postinhibitory rebound (PIR), for the integrators that are typically used. This results in much greater robustness to heterogeneity that actually increases as the average participation in spikes per cycle approximates physiological levels. Moreover, dynamic clamp experiments that show autapse-induced firing in entorhinal cortical interneurons support the idea that PIR can serve as a network gamma mechanism. Furthermore, parvalbumin-positive (PV^{+}) cells were much more likely to display both PIR and autapse-induced firing than GAD2^{+} cells, supporting the view that PV^{+} fast-firing basket cells are more likely to exhibit class 2 excitability than other types of inhibitory interneurons.

**SIGNIFICANCE STATEMENT** Gamma oscillations are believed to play a critical role in information processing, encoding, and retrieval. Networks of inhibitory interneurons are thought to be essential for these oscillations. We show that one class of interneurons with an abrupt onset of firing at a threshold frequency may allow more robust synchronization in the presence of noise and heterogeneity. The mechanism for this robustness depends on the intrinsic resonance at this threshold frequency. Moreover, we show experimentally the feasibility of the proposed mechanism and suggest a way to distinguish between this mechanism and another proposed mechanism: that of a stochastic population oscillator independent of the dynamics of individual neurons.

## Introduction

Several important cognitive functions (Fries, 2009; Lisman and Jensen, 2013) have been hypothesized for gamma oscillations in the 30–90 Hz frequency band (Bressler and Freeman, 1980; Buzsáki and Wang, 2012). Multiple gamma mechanisms have been identified *in vitro* (Bartos et al., 2007; Whittington et al., 2011), including pyramidal-interneuronal network gamma (PING), which requires interplay between populations of excitatory and inhibitory neurons, and interneuronal network gamma (ING), in which purely inhibitory interactions are sufficient. Multiple gamma mechanisms are also likely operative *in vivo* (Wang, 2010; Buzsáki and Wang, 2012). There are two well established theoretical mechanisms that can generate inhibitory network oscillations. In one, the individual neurons are mean-driven, spontaneously spiking neurons in which noise does not affect the mean rate; in the other, individual neurons are fluctuation-driven, Poisson-process-like neurons (Schreiber et al., 2009) driven by noise to fire at random. In the mean-driven regime, individual neurons act as pacemakers that either synchronize or desynchronize their activity depending upon how these oscillators interact when coupled (Wang and Buzsáki, 1996). In the fluctuation-driven regime, individual neurons fire irregularly, but the population rate oscillates (Brunel and Wang, 2003) due to strong, delayed, recurrent inhibition. In the latter case, the network is a stochastic population oscillator. Bartos et al. (2007) argued that the stochastic population oscillator model alone cannot account for the high participation rates of interneurons, which fire on most cycles during both ING (Whittington et al., 2000) and PING (Hájos et al., 2004), and that synchrony in the coupled oscillator network model is destroyed when physiological levels of heterogeneity in excitatory drive are included (Wang and Buzsáki, 1996).

Model neurons used previously to implement the neural oscillator or population oscillator models shared the same excitability type (Hodgkin, 1948; Ermentrout, 1996; Izhikevich, 2007): integrators with a frequency/current relationship allowing arbitrarily slow spiking near threshold. There is another excitability type in which neurons are instead resonators with a cutoff frequency below which they cannot fire. We call this type of resonance “spiking resonance” (Beatty et al., 2015) to avoid confusion with the “subthreshold resonance” (Hutcheon and Yarom, 2000) that is often measured using a ZAP signal. There is strong evidence that the parvalbumin-containing, fast-spiking (FS) basket cells are critical for gamma rhythmogenesis both *in vitro* (Gulyás et al., 2010) and *in vivo* (Cardin et al., 2009). There is also strong evidence that cortical FS interneurons are resonators. Mancilla et al. (2007) found that these neurons rarely fire below 25–30 Hz and Erisir et al. (1999) reported a sharp threshold for the onset of repetitive spiking that is characteristic of a resonator. Tateno et al. (2004) showed several examples of neurons that randomly alternate between regular firing and subthreshold oscillations near the threshold; this “stuttering” is clearly characteristic of a resonator (Izhikevich, 2007). FS hippocampal basket cells have not been studied as carefully, but an earlier study (Pike et al., 2000) revealed that they spiked preferentially in response to sinusoidally varying inputs at 30–50 Hz. Moreover, resonator characteristics and spiking resonance at 40 Hz has been reported for striatal FS cells (Sciamanna and Wilson, 2011).

Recent computational studies (Baroni et al., 2014; Moca et al., 2014) found that resonant interneurons enhance the ability of both ING and PING networks to maintain a uniform frequency in the presence of variable levels of extrinsic drive. Here, we focus on minimal, purely inhibitory ING models and demonstrate robust synchronization in the presence of noise and heterogeneity in synaptic connectivity, resonant frequency, and conduction delays regardless of whether individual neurons are biased in the mean or fluctuation-driven regimes. Furthermore, we explain that the basis for this robustness is spike subtraction resulting in sparse firing and use the dynamic clamp (Economo et al., 2010; Lin et al., 2010) to provide proof of principle for the contribution of resonance to network synchronization.

## Materials and Methods

##### Simulations of resonator neurons.

The two equations that describe the state of each canonical spiking model neuron (Izhikevich, 2003) are as follows:
where the index *j* indicates the *j*th neuron in a population of 300 neurons, time is in milliseconds, membrane potential *v*_{j} is in millivolts, the applied current *I*_{j} is given in nanoamperes and is set to 0.15 unless otherwise noted, and the recovery variable *u*_{j} is dimensionless. We have added a scale factor *k*_{j} to some simulations to take into account that, in addition to variability in external drive, represented by *I*_{j}, the intrinsic dynamics will also be heterogeneous across the population, meaning that, in a biological network, the *F*/*I* curve of all resonant neurons is not identical. The parameters of the resonator were set to *a* = 0.1 and *b* = 0.26, obtained from the MATLAB code at www.izhikevich.org for a resonator cited in Figure 1*K* of Izhikevich, 2004, and Figure 8*K* of Izhikevich, 2007. *v*_{j} is reset when the potential reaches threshold and the recovery variable *u*_{j} exhibits spike-triggered adaptation with *d* = −1 as follows:
The reset potential *c* was changed from −60 mV in the original parameter set to −65 mV because resetting to −60 mV killed the oscillations in the bistable regime described in the Results.

The variable external drive *I*_{j} represents an average general level of excitation and remains inside the brackets for intrinsic dynamics to keep the bifurcation structure uniform. In the network, we have added a noise term and synaptic connectivity outside of the brackets to separate these effects, so Equation 2 is unchanged but Equation 1 becomes Equation 4 as follows:
where *J*_{j}(*t*) is a random process with zero mean and variance σ_{N}, the parameters for the synaptic connectivity term are *E*_{syn} for the synaptic reversal potential, *g*_{ij} for the synaptic conductance, *s*_{ij}(*t*) for the synaptic activation levels, and Δ_{ij} for an optional conduction delay. The noise process was only sampled every 0.1 ms and linearly interpolated between these times to produce consistent results regardless of the time step. *E*_{syn} is set to −70 mV to simulate inhibition and the value of *g*_{ij} in the network is only nonzero for 40 randomly chosen inputs from interneuron *i* to each interneuron *j* for sparsely connected networks (see Figs. 4, 5, 7, and 8) or for all 299 connections for all-to-all connected networks (see Figs. 9, 10). Note that, for both connectivity patterns, autosynapses are prohibited; that is, *g*_{ij} = 0 if *i* = *j*. Additional differential equations are required to describe activation of the biexponential synapses as follows:
The values for the rise and fall time constants were 2 and 5 ms, respectively, unless otherwise noted. The parameter *f* is always set to the value that causes the peak of the biexponential synaptic activation to be 1. *g*_{ij} is set to 0.03 nS for the resonator unless otherwise noted. Autapses were used only in the single neuron simulations in Figures 2 and 3. For the measurement of the phase-resetting curve, the average of the five cycle periods before the untriggered input was used to estimate the intrinsic period. The phase resetting was plotted as the network period divided by the estimated intrinsic period such that delays produced positive resetting.

##### Simulations of integrator neurons.

We model integrator neurons as type 1, conductance-based Wang–Buzsáki model neurons (Wang and Buzsáki, 1996). The current balance equation for each Wang–Buzsáki model neuron is as follows:
Where the capacitance *C* = 1 μF/cm^{2}, *V* is the cell membrane voltage in millivolts and *t* is time in milliseconds. The leak current is given by *I*_{L} = *g*_{L}(*V* − *EL*). The sodium current is given by *I*_{Na} = *g*_{Na}*m*_{∞}^{3}*h*(*V* − *E*_{Na}). The steady-state activation *m*_{∞} = α_{m}/(α_{m} + β_{m}) where α_{m}(*V*) = −0.1 (V + 35)/{exp[−0.1(V + 35)]−1} and β_{m}(*V*) = 4 exp[− (V + 60)/18]. The rate equation for the inactivation variable *h* in the expression for sodium current is as follows: *dh*/*dt* = Φ{α_{h}(*V*)(1 − h) − β_{h}(*V*)*h*} where Φ = 5. The rate constants for the inactivation variable h are given by α_{h}(V) = 0.07 exp[− (*V* + 58)/20] and β_{h}(*V*) = 1/{exp[−0.1(V + 28)] + 1}. The potassium current is given by *I*_{K} = *g*_{K}*n*^{4}(*V* − *EK*), where the activation variable *n* satisfies the following equation: *dn*/*dt* = Φ{α_{n}(*V*)(1 − *h*) − β_{n}(*V*)*n*} where the rate constants for *n* are α_{n}(*V*) = −0.01(*V* + 34)/{exp[−0.1(V + 34)]−1} and β_{n}(*V*) = 0.125 exp[−(*V* + 44)/80]. The reversal potentials *E*_{Na}, *E*_{K}, and *E*_{L} were set to 55, −90, and −65 mV, respectively. The maximal sodium (*g*_{Na}), potassium (*g*_{K}), and leak (*g*_{L}) conductances were set to 35, 9, and 0.1 mS/cm^{2}, respectively. *I*_{i} is the applied current that functions as the external drive. Unless otherwise stated, the values for the various parameters were equal to those given above. The synaptic current is as described for the Izhikevich model except that *E*_{syn} is equal to −75 mV. A synaptic event is triggered by an upward crossing of a 10 mV threshold.

##### Network simulations.

Both models were implemented in the simulation package NEURON (Hines and Carnevale, 1997). For systems with negligible delays, we often set the delays to 0.1 ms for greater computational efficiency while using multithreaded/MPI parallelization. Synaptic activation initially was set to zero for all simulations.

In Figures 4, 5, 7, and 8, state variables *v* and *u* were randomly initialized from Gaussian probability distributions with −51.86 ± 20 mV and −15 ± 5 (mean ± SD), respectively.

In Figure 8*C*, individual synaptic conductances were picked from log-normal distribution (Song et al., 2005). The variance of the distribution was increased by increasing the scale parameter of the distribution while also adjusting the location parameter to keep the mean of the distribution equal to 0.03 nS. In Figures 9 and 10, initial conditions were *v* = −65 mV and *u* = −16.5 for each resonator and *v* = −15 mV, *n* = 0, and h = 1 for each integrator to ensure that each neuron generated at least one spike. Moreover, in Figures 9 and 10, we used all-to-all connectivity instead of 40 connections per neuron to remove spatial inhomogeneity from the network; this manipulation favors the integrators because the resonators are more robust to spatial inhomogeneity and therefore perform even better than the figures indicate.

In Figure 9, a two-step procedure was used to set up a fair comparison between resonators and the integrators in the tonic regime. First, the external drive parameter *I*_{j} was set to give a 32–36 Hz spontaneous firing rate in isolated neurons of both model types (i.e., 5 μA/cm^{2} for integrator and 0.2 nA for resonator). Then the synaptic conductance *g*_{syn} was adjusted so that both networks had a frequency of 25–27 Hz (i.e., 1 nS/cm^{2} for integrator and 0.035 nS for resonator). A 1 ms conduction delay was required to get global synchrony in the integrator network, so we added a 1 ms delay to both networks. A different strategy was used to devise a fair comparison between the resonator network and the integrator network in the phasic regime. Because the synaptic conductance should be relatively strong in the phasic regime, the excitatory/inhibitory balance was preserved across networks with *g*_{syn}(*v*_{rest} − *E*_{syn}) = 0.125 *I*_{j} for both models (i.e., *g*_{syn} = 0.125 μS/cm^{2} for integrator and 0.005 nS for resonator). In both cases, we then varied *I*_{j} across the population of integrators by picking each value from a Gaussian distribution. We scaled the intrinsic dynamics across the population of resonators by varying *k*_{j} to get approximately the same distribution of intrinsic frequencies. The coefficient of variation (CV) for intrinsic frequencies was then calculated from mean and SD in current using *F–I* curves for both model types.

In Figure 10*A*, *I*_{j} was set below threshold to 0.15 nA for resonator and 1.5 μA/cm^{2} for integrator. Noise amplitude was set to give a mean firing rate of 10 Hz in both isolated model neurons (i.e., 3 μA/cm^{2} for integrator and 0.87 nA for resonator). Then, *g*_{syn} was set to produce the same level of synchronization for both networks (0.07 μS/cm^{2} for integrator and 0.03 nS for resonator). We then varied the total synaptic conductance across the population by using these values as the mean and picking the value of *g*_{syn} for all synapses onto a given neuron from a Gaussian distribution truncated to eliminate negative values (Fig. 10*A*). For Figure 10*B*, *I*_{j} and *g*_{syn} had the same values as in Figure 10*A*, with *g*_{syn} set everywhere to its mean value, but the noise amplitude was swept through the range [2.8, 90] μA/cm^{2} for integrator and [0.06, 13] nA for resonator. We have summarized noise and connectivity parameters for all figures in Table 1.

##### Novel synchrony measure.

Oscillatory frequency can vary over time in a nonstationary fashion (Siapas et al., 2005; Thounaojam et al., 2014), which could cause linear methods such as cross-correlograms to underestimate synchrony. The vector strength statistic or its square, the *R*^{2} statistic from circular statistics, is often defined in terms of the average period (Batschelet, 1981), which could also underestimate synchrony. We have chosen to quantify synchrony used a method we developed based on the cycle by cycle population period, as defined by the peaks of the population histogram. The spikes were binned in 1 ms windows and a pulse with the height determined by the number of spikes in the bin was placed at the center of the bin. The resultant pulse train was low-pass filtered by convolution with a Gaussian kernel with σ = 10 ms and a length of 100 ms, which produced a time series with clear peaks in the network activity to use as a clock to compute the vector strength. The peaks were detected using downward zero crossings of the first derivative of the smoothed signal. The peak at the beginning of each network cycle was assigned a phase value of 0 and, at the end, a phase value of 2π. Each spike was assigned a phase depending upon where it fell within the cycle and a vector length of 1 starting from the origin. Then the vector sum of all the spike vectors constructed in this fashion was taken and normalized by the number of vectors, resulting in a vector with an average phase and a length between 0 and 1. The *R*^{2} statistic is the square of the vector length. This measure quantifies the level of synchrony of the individual neurons with the population rhythm, but gives no information regarding the average level of participation in the oscillation, which we track separately as the average number of spikes per cycle (SPC) normalized by the population size. Note that random peaks of spike rate can be detected in finite networks even in the absence of network oscillations, for example, in random firing or in phase dispersion. In this case, the *R*^{2} vector length is close to 0, but the SPC still reflects the mean firing rate of the population. All simulations were 10 s in duration. Simulations were classified as nonoscillatory if the average spike time was <2.5 s, indicating that spiking dies out, and were rejected if the average spike time was >7.5 s, indicating an unacceptably long period of transient activity. We doubled simulation times and cutoff values and obtained similar simulation statistics, verifying the effectiveness of these criteria.

##### Experimental methods.

All electrophysiological experiments were conducted according to protocols approved by the University of Utah Animal Care and Use Committee. Brain slices were harvested from 18- to 35-d-old *cre*-dependent GAD2-IRES-tdTomato (Taniguchi et al., 2011) or PV-tdTomato (Hippenmeyer et al., 2005) transgenic mice, which labeled all glutamic acid decarboxylase 2 gene (GAD2)- and parvalbumin (PV)-expressing cells, respectively. All interneurons recorded were located in layers 1, 2, and 3 of dorsal medial entorhinal cortex. All PV^{+} cells were in layers 2 and 3. These mice were anesthetized with isofluorane and decapitated. The brain was then harvested, chilled in sucrose-substituted artificial CSF (ACSF) containing the following (in mm): 185 sucrose, 2.5 KCl, 1.25 NaH_{2}PO_{4}, 10 MgCl_{2}, 25 NaHCO_{3}, 12.5 glucose, 0.5 CaCl_{2}), and cut parasagittally into 300-μm-thick slices using a vibrating microtome (Vibratome VT1200; Leica). Slices were incubated for 15 min in ACSF containing the following (in mm): 125 NaCl, 2.5 KCl, 1.25 NaH_{2}PO_{4}, 10 MgCl_{2}, 25 NaHCO_{3}, 25 glucose, 2 CaCl_{2}) at 37°C, and then allowed to recover for at least 30 min at room temperature. For recordings, slices were transferred to a heated (32–34°C) slice chamber (Warner Instruments) mounted on an upright microscope stage (BX53; Olympus) and perfused with 95/5% O_{2}/CO_{2} ACSF. GAD2/tdTomato^{+} and PV-tdTomato^{+} neurons were visualized using fluorescence and whole-cell patch clamped using patch pipettes (5–6 MΩ) fabricated from borosilicate glass (1.5 outer diameter, 1.1 inner diameter; Sutter Instruments) and filled with artificial intracellular fluid (ICF) containing the following (in mm): 120 K-gluconate, 5 MgCl_{2}, 0.2 EGTA, 10 HEPES, 20 KCl, 7 di(tris) phosphocreatine, 4 Na_{2}ATP, 0.3 Tris-GTP. Presented data were not corrected for the junction potential, which was assumed to be 10–12 mV.

For all experimental protocols, synaptic conductances were simulated using dynamic clamp software (Dorval et al., 2001; Bettencourt et al., 2008; Lin et al., 2010) on a Pentium 4 computer running Linux Ubuntu with a patched version of the real-time application interface (RTAI) kernel. Voltage was measured and a control current applied with a MultiClamp 700B amplifier (Molecular Devices
where *g*_{max} is the maximal conductance, *V* is membrane voltage, *E*_{syn} is the reversal potential of the synapse (−75 mV for inhibitory, 0 mV for excitatory), and *s*(*t*) is the difference of two exponentials with time constants of τ_{rise} = 1 ms and τ_{fall} = 3 ms. Maximal conductances ranged between 0.8 and 2 nS. To examine a neuron's ability to exhibit post-inhibitory rebound firing, single artificial inhibitory postsynaptic conductances (IPSGs) were elicited near threshold with varying maximal conductances. If the initial trial elicited little to no postinhibitory rebound (PIR) firing, then the cell was depolarized further and the trial was conducted again. This process was repeated until the cell exhibited faster than 1 Hz firing. The most depolarized IPSG trial was then used to test for PIR firing. This process ensured that PIR firing was not undetected due to the cell being too hyperpolarized. In the most depolarized trial, if a neuron was able to exhibit PIR firing in >15% of induced IPSGs, then it was determined to be capable of PIR spiking.

In neurons exhibiting PIR firing, an artificial autapse condition was simulated wherein the detection of an action potential (determined by a crossing of the −20 mV potential threshold) was followed by an artificial IPSG with a 2 ms delay. Some neurons were able to exhibit IPSG-induced PIR firing, but were not able to maintain autapse-induced firing. At least 10 autapse trials at maximum IPSG amplitude and varying membrane potentials were conducted to test for autapse-induced firing. This process ensured that small IPSG inputs or holding membrane potentials did not produce false negatives in determining whether cells could exhibit autapse-induced firing. To show phase resetting by an additional non-autapse-induced IPSG, in some trials, an IPSG was inserted randomly 100–300 ms after the initiation of autapse-driven firing. The magnitude of the disruptive IPSG was matched to an initial, PIR-inducing IPSG and was generally 100–150% of the magnitude of the IPSGs used to maintain the autapse. This ensured that the disruptive IPSG was large enough to trigger an IPSG. For the phase-resetting curve, the median interspike interval (ISI) for all spikes in the autapse-driven regime was used to estimate the intrinsic period. The phase resetting was plotted as the total ISI of the cycle that received the disruptive IPSG divided by the estimated intrinsic period, such that delays produce positive resetting. A one-tailed *Z* test was used to determine whether the assumed binomial distributions for exhibiting PIR and exhibiting autapse-induced oscillations had a higher probability for the PV^{+} population compared with the GAD2^{+} population.

In addition, the likelihood of the presence of a discontinuity in the input current-to-firing frequency relationship was assessed by measuring this relationship in current clamp. In these trials, a bias current was injected to hyperpolarize the cell to −70 and a series of 1 s current pulses were injected between −100 pA up to peak firing potential (between 500 pA and 1500 pA, depending on the cell) in 20 pA increments. The lowest non-zero firing frequency was reported using this protocol, with a higher value corresponding to greater likelihood of a discontinuity associated with type 2 excitability.

## Results

### Resonant interneurons

Postinhibitory rebound can occur if an inward current is activated or deinactivated by hyperpolarization or, conversely, if an outward current is deactivated or inactivated by hyperpolarization. The net, short-term effect of the hyperpolarization is to shift the balance in favor of the inward currents if the original membrane potential is revisited before the recently engage PIR current returns to steady state. We are interested here only in PIR mechanisms fast enough that can be engaged by inhibitory postsynaptic potentials (IPSPs) generated at GABA_{A} synapses; for example, removal of sodium channel inactivation. In a type 2 oscillator, the net steady-state current at membrane potentials traversed during the ISI is outward, so pacing can only be maintained when the ISI is traversed quickly enough that the inward current evoked by the after-hyperpolarization can persist (Prescott et al., 2008); this rate sets the lower bound on spiking frequency. Otherwise, the oscillation will stop at the point the inward and outward currents come into balance. Therefore, in a sense, pacemaking in type 2 neurons is driven by PIR. One signature of a type 2 neuron is an oscillatory tendency near spike threshold. The membrane potential of a quiescent neuron generally relaxed exponentially back to rest after a perturbation (Kandel et al., 2000). Conversely, a type 2 neuron can exhibit damped oscillations that overshoot the holding potential after the perturbation: this overshoot is evidence of rebound. In the model, the slow adaptation variable is reduced by hyperpolarization, which favors the balance of inward currents so a simulated IPSP produces a damped oscillation (Fig. 1*A1*). Moreover, an IPSP of sufficient amplitude produces a rebound spike (Fig. 1*A2*).

To test whether GABAergic interneurons could exhibit type 2 dynamics and postinhibitory rebound, we used Linux-based dynamic clamp (Economo et al., 2010; Lin et al., 2010) to apply simulated virtual conductance waveforms that would be expected during IPSPs mediated by GABA_{A} receptors. Measurements were made in GAD2^{+} neurons, presumably representing the entire GABAergic population, and in PV^{+} neurons, representing mainly FS basket cells that play an important role in gamma. A key observation is that, if the neuron is too hyperpolarized, the PIR mechanism may be already fully engaged so additional hyperpolarization has no effect. Therefore, no rebound was observed from rest. However, when the membrane was slightly depolarized by constant current injection, an overshoot of the membrane potential accompanied by a damped oscillation could be indeed be detected (Fig. 1*B1*) in many neurons. When the amplitude of the simulated conductance was increased, a rebound action potential was triggered (Fig. 1*B2*) for 15% or more of trials in 31 of 81 (39%) of GAD2^{+} cells and 22 of 34 (65%) of PV^{+} cells. This effect is significant (one-tailed *Z* test, *p* < 0.01). Keeping in mind that the PV^{+} population is a subset of the GAD2^{+} population, these results indicate that PIR-induced firing is very common among FS basket cells, but uncommon among other cell types. Therefore, it is likely that a large fraction of the basket cells in dorsal MEC exhibit type 2 dynamics. Moreover, PV^{+} cells had a much faster lowest possible firing frequency (78.9 ± 6.2 Hz) than GAD2^{+} cells (16.0 ± 3.2 Hz, *p* < 0.001). This suggests that PV^{+} cells had more easily observable *F–I* curve discontinuities, likely associated with type 2 excitability, compared with GAD2^{+} cells.

### Proof of principle PIR network dynamics

We next wanted to test whether a constant-latency PIR could form the basis for a network oscillation. We first used a single, isolated model neuron and elicited a PIR spike with an initial input (triangle in Fig. 2*A1*). Each spike in the model neuron triggered a synaptic conductance waveform. We reasoned that this self-feedback signal could approximate the population spike evoked in other interneurons by network feedback when the model neuron was embedded in a network. This protocol produced an oscillation sustained by feedback (Fig. 2*A1*). Current noise was added to the model neuron to give a spread in the ISI histogram (Fig. 2*B1*). We then repeated this protocol in a biological setting, again using the virtual GABA_{A} conductance implemented using the dynamic clamp. In this case, after an initial IPSP to evoke PIR, subsequent feedback IPSPs were triggered by each spike emitted by the interneuron and, again, a sustained oscillation was observed (Fig. 2*A2*). Indeed, the histogram of ISIs observed in the self-connected neuron was narrow (Fig. 2*B2*), indicative of a fairly constant latency, with a peak in the low gamma range. Overall, we observed autapse-supported periodic gamma in only 9 of 81 (11%) of GAD2^{+} cells, but in 15 of 34 (44%) of PV^{+} cells. This difference between PV^{+} cells and the broader GAD2^{+} population was highly significant (one-tailed *Z* test, *p* < 0.001), consistent with the interpretation that only the FS interneurons are able to support ING via resonance and PIR-induced firing. Notably, the PV^{+} cell population that exhibited rebound action potentials had a significantly shorter latency (input-to-spike delay) than the GAD2^{+} cell rebound-firing population (25.4 ± 1.8 ms vs 45.5 ± 2.9 ms; *p* < 0.001). PV^{+} cells that could also exhibit the autapse-induced firing had significantly shorter latencies than those that could not (22.2 ± 1.7 ms, *n* = 15, vs 31.6 ± 2.9 ms, *n* = 7; *p* < 0.01). A longer latency allows more opportunity for noise to kill the oscillation.

Using a model neuron with an autapse as in Figure 2*A*, an additional, untriggered input was applied after a feedback IPSP, but before the corresponding rebound spike was emitted (triangle in Fig. 3*A*). The phase response curve (PRC) for the autapse-induced oscillation was measured by defining a single cycle as the ISI, then measuring the phase delay resulting from the additional input at each phase in the cycle. We hypothesized that the untriggered input would erase the memory of the most recent autapse inhibition and effectively restart the clock for the PIR latency between the hyperpolarizing input and the next rebound spike. If the input is assumed to be saturating, then delivering an additional input at approximately the same time as the autapse should have no effect and, indeed, we observed very little phase resetting at early phases (Fig. 3*A2*). Moreover, if the latency were truly constant, then the perturbed cycle period is simply the sum of the time elapsed between the previous spike and the input and the constant latency equal to the period of the autapse-induced oscillation. Therefore, the phase delay is predicted to be exactly equal to the phase and the PRC should be a line passing through the origin with a slope of 1. Indeed, we obtained a PRC with a slope of 1.1 in the model. In the model, sometimes an input applied during the upstroke of the action potential could not prevent the action potential from occurring, resulting in a branch with little to no resetting at late phases. This branch was ignored in the slope calculation.

An additional, untriggered input was applied (triangle in Fig. 3*B*) to measure the PRC for autapse-induced oscillations in 10 PV^{+} cells (Fig. 3). We obtained a linear relationship of slope 1.32 ± 0.12 with a coefficient of determination of 0.83. Therefore, the observed latency is close to the theoretical prediction in terms of the shape, but is not exactly constant. As we will see in the next sections, the fact that the slope is greater than required for a constant latency does not substantially interfere with the robustness of the network oscillation. What is important is that phase-resetting data points do no occur frequently very far below the diagonal with a slope of 1 because these points with short latencies would lead to spikes before the rest of the population. The next network wave of inhibition will reset the phase of neurons on track to produce a spike with a longer latency than the network period, so points far above the diagonal are not problematic (see next sections). Moreover, our simulations (data not shown) revealed that, for the minimal input required to evoke PIR, the slope was larger than 1 (up to 1.3), but as the strength of the synapse was increased, the slope converged to 1.

### Cycle skipping in a noisy network

Having established the physiological feasibility of an ING mechanism specific to type 2 neurons, we proceeded to use the model to demonstrate the high levels of robustness to heterogeneity that can be achieved by this network mechanism. First, we introduced a type of heterogeneity in the time domain by applying independent current noise to each model neuron (Fig. 4). In the absence of noise, global synchrony is an emergent property of this network that arises from symmetry (Golubitsky and Stewart, 2006). In this particular case, the individual neurons are biased in the subthreshold (nonoscillatory) regime, but random initialization causes random neurons to emit spikes, recruiting the population into a globally attracting PIR based network oscillation (for *g*_{syn} = 0.03 nS). In the absence of any noise or heterogeneity, all neurons fire at the exact same time because they are identical and receive an identical number of identical inputs on each network cycle. No neuron is prevented from firing because inhibition is not instantaneous, but in general has at least a small conduction delay, and a finite rise time is required for inhibition to reach full strength.

When uncorrelated noise is applied to each neuron in the network (bar in Fig. 4), jitter emerges in the firing times of the interneurons during a population burst. Neurons that, in the absence of coupling, would have spiked toward the end of the population burst instead have their latency reset as in Figure 3 and do not fire on that cycle (Fig. 5*A1*). The trace of the membrane potential in a representative neuron (Fig. 4*A1*), as well as the raster plot (Fig. 4*A2*), show that cycles are skipped. The spike time histogram (Fig. 4*A3*) shows that the firing rate drops dramatically in the presence of noise, but cycle skipping allows synchrony to be fairly well preserved because mistimed spikes are largely precluded. Finally, the histogram of all ISIs (Fig. 4*B*) for each neuron within the network includes a peak at the network frequencies but also peaks at integral multiples, with integral multiplier indicating the number of cycles that were skipped by individual neurons. The connectivity in this network is not all-to-all, but instead is sparse and random. Although, in this simulation, all neurons have the same number of postsynaptic synapses, a random number of these inputs are active due to the participation of random neurons in each cycle. This is a critically important consequence of spatial inhomogeneity in the subsets of presynaptic neurons that is unmasked by the noise. In this particular simulation, as stated above, all neurons are biased in the subthreshold regime; therefore, in this case, randomness is increased because some neurons on each cycle may not receive sufficient hyperpolarizing input to evoke a rebound spike on the next cycle.

### What controls the window in which firing is allowed?

To explain how cycle skipping enforces sparse synchrony, we compared simulations in coupled versus uncoupled networks (Fig. 5*A*). The networks were randomly initialized as described in the Materials and Methods, and then a strong common inhibitory input was applied. The white region of the peristimulus time histogram shows the dispersion of the spikes in an uncoupled, sparsely connected but otherwise homogeneous, noisy network. The gray region shows the dispersion of the spikes in a coupled network: the gray region is a subset of the white region in which spikes that occurred later were preferentially skipped. Increasing the conduction delay from 0 to 3 ms broadened the window in which spikes could be fired (Fig. 5*B*) by delaying the arrival of the barrage of inhibition.

### Dynamical systems explanation of cycle skipping that underlies sparse synchrony

What determines whether a particular neuron skips a cycle? The reason that the neurons either spike with the population or wait until the next cycle to have another chance to fire results from the qualitative nature of the dynamics near a subcritical Hopf bifurcation, which we will explain pictorially by potting the two state variables, *v* and *u*, in the (*v*, *u*) plane. For a subthreshold, quiescent neuron (Fig. 6*A*), there are two points at which the net current is 0 and the recovery variable *u* is at its steady-state value: the stable resting potential of the model neuron (open circle) and a second point (filled diamond) at which these conditions are also met. The stable resting point is a focus, meaning that trajectories spiral into it. The spiral in the phase plane produces the damped oscillations seen in Figure 1, *A1* and *B1*. Conversely, trajectories in the phase plane are attracted to the saddle point (without spiraling) from one direction (the stable manifold) and repelled in another (the unstable manifold). The repulsion results from the positive feedback built into the model by positive sign of the coefficients of the *v* and *v*^{2} terms in Equation 1, so any tiny deviation from the saddle point results in the upstroke of an action potential. The PIR dynamics can be simply explained by stating that inhibition pushes the trajectories toward the left. After a reset, trajectories move toward the open circle. If inhibition is not sufficient to push the trajectory outside of the boundary to the left, then no action potential is generated and a cycle is skipped. Conversely, if inhibition pushes the trajectory to the left outside the boundary, then it lands in a confined region of the phase plane from which there is an approximately constant latency until an action potential is generated unless inhibition from the next population burst causes it to cross the boundary before a spike is generated and skip a cycle.

For values of the applied current between 0.1795 and 0.2625 nA, an individual resonator model neuron is bistable between quiescence and spontaneous pacemaking; this is characteristic of the dynamics near a subcritical Hopf bifurcation. As the applied current is increased from values that produce quiescence, the vector field “pinches off” to form a closed curve that is an unstable limit cycle (dashed curve in Fig. 6*B*); this limit cycle now assumes part of the role of marking the border between generating a spike or skipping the cycle. Trajectories that remain inside or reenter the limit cycle spiral into the stable rest potential and skip the next cycle. However, the remainder of the boundary is formed by the membrane potential nullcline (gray curve), which divides the plane into regions in which the membrane potential is decreasing (leftward movement) or increasing (rightward movement). Trajectories that fall above the limit cycle but to the left of the right branch of the nullcline also skip a cycle. During the decay of inhibition, the trajectories begin to move to the right, and trajectories to the right of the nullcline or below the limit cycle when the effect of inhibition has dissipated fire an action potential. In general, neurons that fire early in the cycle avoid being pushed across the boundary that causes a missed cycle.

In the suprathreshold regime (Fig. 6*C*), the rest potential becomes unstable in a way that forces trajectories to spiral away from it, such that model neurons become spontaneous pacemakers. Moreover, the stable manifold of the saddle point now starts at the now unstable “resting potential” (filled circle) and continues through the saddle point via the unstable manifold to infinity. Trajectories to the right of below this boundary result in an action potential, whereas others result in a skipped cycle. Each of these three scenarios (subthreshold, bistable, suprathreshold) indicates an oversimplified static picture of the dynamics because the level of inhibitory conductance experienced by each neuron is constantly changing, resulting in a movie rather than a snapshot (Rotstein, 2015). However, in all cases, one of these snapshots qualitatively characterizes the dynamics; therefore, the mechanism of cycle skipping is robust and conceptually similar in all regimes. Although this model has simplified dynamics, it is dynamically equivalent to the more complex and better known FitzHugh–Nagumo and Hodgkin–Huxley models and the results presented in the networks in this study should be general.

### How does inhibition control the spiking window?

We next examined the effect of the synaptic parameters on the width of the window in which firing was allowed (Fig. 7*A*). We assessed by width of the window using the SPC. We fixed τ_{rise} at 2 ms and allowed τ_{fall} to increase from 2 to 10 ms and allowed *g*_{syn} to vary from 0 to a very strong value. In this parameter space, there is a region of silence, a region of sparse synchrony characterized by *R*^{2} > 0.9 and a region of intermittent synchrony in which *R*^{2} drops off rapidly. The latter two regions result from random initialization: a globally synchronous solution is bistable with sparse synchrony. Because there is no noise in this figure, the only source of variability is the sparse connectivity. Spatial inhomogeneity is sometimes called quenched randomness (van Vreeswijk and Sompolinsky, 1998). With all to all connectivity, sparse synchrony (SPC <1) is not exhibited in the absence of independent noise sources or other forms of heterogeneity. With full connectivity and no variability, global oscillatory synchrony (SPC = 1) or the trivial version of synchrony in which all neurons are silent (SPC = 0) are the only solutions exhibited by the network. With sparse connectivity, one or both of these solutions can be bistable with sparse synchrony. However, in Figure 7, we are only interested in characterizing how the synaptic parameters affect width of the window in which firing is allowed, so we ignore this bistability.

For *g*_{syn}, there a minimum value below which PIR is not generated (Fig. 1) but it is barely visible on the scale in this figure. As *g*_{syn} is increased, the window for spiking becomes narrowed, as evidenced by the decrease in SPC and by the appearance of multiple peaks in the ISI histograms (Fig. 7*B*). Peaks other than the first one indicate skipping one or more cycles and broader peaks reflect a broadening of the window. Increasing *g*_{syn} increases the degree of quenched randomness and destroys the global attractiveness of full synchrony observed in the noiseless part of Figure 4*A*. For an instantaneous conductance (delta function), there is no upper bound on *g*_{syn} because the effect of the synaptic conductance saturates as the membrane potential approaches the reversal potential. However, synchrony is weakened as *g*_{syn} is increased because the inhibition becomes essentially tonic or, on much of the time, blurring the simple picture presented in Figure 6 in which a neuron either skips a cycle or it does not. The simulations are not stationary in the hatched regime in Figure 7*A*, instead sparse synchrony waxes and wanes, with less synchrony as *g*_{syn} continues to increase. This effect of increasing *g*_{syn} is, of course, exacerbated as τ_{fall} is increased. τ_{fall} has no lower bound because it is not allowed to drop below τ_{rise}, but it clearly has an upper bound when inhibition is so long-lasting that the stable fixed point in Figure 6*A* is always achieved rather than PIR.

We did not systematically investigate the effect of τ_{rise} because increasing τ_{rise} has a similar effect as increasing the conduction delays in that it broadens the window. However, we did perform simulations with single exponential synapses that were allowed to reach their maximum value with no delay after the presynaptic neuron reached the spike threshold. For a fully connected network, very fast inhibition disrupted synchrony, in agreement with decades of literature (Van Vreeswijk et al., 1994). Global synchrony was also disrupted in a sparsely connected network wired as described in the Materials and Methods. However, sparse synchrony persisted due to cycle skipping. Global synchrony requires that some time elapse before a spike in one neuron delays the occurrence of a spike in its neighbors, otherwise global synchrony is destabilized. Cycle skipping abolishes this requirement and sparse synchrony persists in the face of strong, instantaneous inhibition.

### Robustness of synchrony in the presence of different types of heterogeneity

We next systematically explored the effect of different types of heterogeneities, starting with current noise (as in Fig. 4) that provides heterogeneity in the temporal domain. As the noise level increased, the SPC (dashed black curve in Fig. 4), which is the fraction of network cycles in which each neuron participated on average, decreased to maintain synchronization in a fashion similar to that observed in Figure 4*A*. The level of the current noise applied to each model neuron was varied from 0 to a level (1.8 nA) that produced a high level of spontaneous firing (∼22 spikes/s) in otherwise quiescent neurons (Fig. 8*A*). The *R*^{2} vector strength (solid black curve) measure of synchronization remained above 0.7 even for this strong level of noise. Next (Fig. 8*B*), we varied the level of external, constant excitatory current, or drive across the network, introducing spatial heterogeneity. The network responded to spatial heterogeneity in exactly the same fashion that it responded to heterogeneity in the time domain (i.e., noise). The vector strength dropped off gradually as the CV of the drive was increased to 1 and the SPC again decreased to allow the network to maintain synchrony by preferentially dropping misaligned spikes. The network spanned the range of oscillators and quiescent neurons, with an average drive of 0.2 nA compared with a level of 0–0.1795 nA for quiescence, 0.1795–0.2625 nA for bistability between quiescence and pacemaking, and pacemaking only at drives >0.2625. We then introduced spatial inhomogeneity in a different way by randomizing the values of individual synaptic conductances in Figure 8*C*, as described in the Materials and Methods and assessing this variability using the CV of the total conductance resulting from the 40 synapses per neuron across the population of 300 neurons. Network synchrony as measured by the *R*^{2} vector strength was exceeding robust to this type of variability, likely because the PIR latency of individual neurons is weakly dependent on the strength of the inhibitory conductance so long as it is sufficient to elicit a rebound spike as shown in Figure 1*B*). Nonetheless, the SPC decreased as heterogeneity in connectivity was increased, but the synchronization remained near perfect for a CV of 1. Finally (Fig. 8*D*), we varied the conduction delays from a mean of 3 ms as the SD of a Gaussian distribution was increased to 14 ms (the distribution was truncated to avoid values of delay <0.1 ms). Once again, synchronization remained high as the variability was increased to very high levels. In contrast to the other types of heterogeneity, for variability in the delays, the SPC only declined slightly until a threshold of a SD of ∼5 ms was reached, then they dropped sharply to allow the network to maintain synchrony. Introducing variability in the conduction delays spreads out the inhibitory barrage in time, even for the case in which all neurons fire simultaneously. For small variability, the inhibitory barrage occurs after the population burst without suppressing any neurons. For an SD of delays >5 ms, the inhibitory barrage was sufficiently spread out that it began to suppress some neurons on each cycle. Therefore, increasing the SD of the delay narrows the window of opportunity to fire.

We also examined how gap junctions affect synchronization and the robustness of the mechanism under study. We inserted an electrical synapse with a conductance of 2 μS in parallel with each chemical synapse in the model. In general, the R2 and SPC values were very similar to those shown in Figure 8 except, in Figure 8*B*, there was a slight increase in *R*^{2} for large SD.

### Comparison of purely inhibitory networks of suprathreshold integrators or resonators

We varied the parameters of each model as described in the Materials and Methods to obtain a distribution in the frequency and plotted the *R*^{2} vector strength as a function of the CV of the frequency. As expected from previous work (Wang and Buzsáki, 1996; White et al., 1998), the ability of a network of integrators to synchronize in the presence of such heterogeneity is very limited. There are two failure modes: phase dispersion, in which *R*^{2} plummeted but participation (SPC) remains high, and suppression, in which *R*^{2} remains high but the faster neurons suppress other neurons so the participation drops. The *R*^{2} vector strength for a network of integrators plummeted near 0 with a CV of only ∼0.1 or 10% (dashed curve in Fig. 9*A*). Conversely, a comparable tuned network of resonators (see Materials and Methods) was able to maintain strong synchronization (solid line near 1 in Fig. 9*A*) at a CV of 0.8, at a cost of decreasing the SPC (Fig. 9*B*). The integrator network was tuned in a phasic regime with stronger inhibition and it failed in a manner similar to the comparably tuned resonator network (see Materials and Methods), but much more precipitously than the resonator network in terms of participation (cf. Fig. 9*B* and *C*2 vs *C3*). However, it is important to note that these simulations were conducted in networks with all-to-all connectivity, which shows the integrators in the best possible light because they do not show the robustness to sparse connectivity of resonator networks. Previous analyses suggested that neurons with type 2 excitability and type 2 PRCs cannot synchronize via inhibitory connections (Hansel et al., 1995; Achuthan and Canavier, 2009). However, we find that strong inhibition that invokes PIR can decouple the shape of the PRC from the excitability type (Ermentrout, 1996) and that this is sufficient to allow global synchronization of neurons with type 2 excitability coupled via inhibition.

### Comparison of purely inhibitory networks of subthreshold integrators or resonators

We then wanted to make a fair comparison between the stochastic population oscillator model and the PIR network or resonators. We took the same network of 300 neurons, made sure that the constant level of drive resulted in a quiescent neuron, and then adjusted the current noise level so that the isolated resonator and integrator model neurons had the same baseline frequency. Then, the synaptic conductance in the network was adjusted so that each connected network had an *R*^{2} vector strength of 0.97. At this value of *R*^{2}, the SPC was ∼0.2 in the resonator network and ∼0.006 in the integrator network. First, we demonstrated that synchronization in the network of resonators was more robust to variability in the connectivity than in networks of integrators, as shown by the more rapid decline of vector strength for the integrator network (dashed line in Figure 10*A*) compared with the resonator network (solid line in Fig. 10*A*) as the CV of the total synaptic conductance was increased. Next, we examined the relationship between participation of individual neurons in the network oscillation as measured by the SPC and synchrony as measured by *R*^{2}. To manipulate the SPC through a range, we varied the noise amplitude as described in the Materials and Methods. In the integrator network, increasing the SPC from 0.006 to ∼0.01 resulted in a precipitous drop in *R*^{2} vector strength (dashed curve in Fig. 10*B*). Conversely, increasing the noise in the network of resonators changed the SPC in the opposite direction, causing a decrease in participation to an SPC of ∼0.1, and also decreased R_{2}, but only slightly, until an inflection point was reached. After the inflection point, *R*^{2} decreased sharply, but participation began to increase. Examination of the ISI histograms suggest that regime in which participation increases with increasing synchronization corresponds to a PIR cycle-skipping oscillation, whereas the regime in which participation increases with decreasing synchronization more closely resembles a stochastic population oscillator. In fact, the stochastic population oscillator composed of integrators shows a similar dependence of decreasing synchrony with increasing participation. In sum, Figure 9 shows that, in the suprathreshold regime, resonators are much better able to synchronize in the face of heterogeneity in frequency than integrators and Figure 10 shows that, in the subthreshold regime, the synchronization of networks of resonators operating under a PIR mechanism is largely insensitive to increasing participation, but that of network of integrators deteriorates rapidly with increasing participation. Only when the noise is large enough to overwhelm the PIR mechanism does the resonator network begin to resemble qualitatively a stochastic population oscillator. We suspect that noise disrupts the vector field, illustrated by the snapshots in Figure 6, which constrains trajectories to either fire together or skip a cycle. Consistent with this idea, we decreased the factor *k*_{j} in Equations 1 and 2 to decrease the strength of the vector field by the same factor at all points and found that the noise required to reach the inflection point in Figure 10*B* decreased approximately proportionately (data not shown).

It is clear in the histograms in Figures 7*B* and 10*C* that, as the robust sparse synchrony due to PIR is lost due to increasing noise or synapse strength, the “forbidden zones” between peaks at multiples of the network ISI fill in. In Figure 6, *A* and *B*, the neurons have a stable fixed point inside of a well defined boundary and, in Figure 6*C*, the situation is somewhat analogous because of the very weak repulsion of the unstable fixed point. Forbidden zones remain as long as the noise is not strong enough to move trajectories from the fixed point across the boundary. To address the effect of increasing synaptic strength in Figure 7, it is important to keep in mind that the boundaries in Figure 6 were drawn for an isolated neuron in the absence of inhibition to explain the trajectories after the inhibition has dissipated. For sufficiently strong inhibition, the inhibition never dissipates, so the phasic effect of inhibition is largely lost and the timing of spikes becomes noise driven.

### Summary

We have demonstrated that interneurons with type 2 excitability are not uncommon in the entorhinal cortex and the literature provides evidence that they can be found in other brain areas as well. We provided a “proof-of-principle” example experiment using the dynamic clamp to show that inhibitory network feedback can produce an approximately periodic rhythm in otherwise quiescent interneurons with type 2 excitability. We showed that inhibitory networks of model neurons with type 2 excitability are remarkably robust to heterogeneity introduced by uncorrelated noise or spatial inhomogeneities introduced in the level of extrinsic drive and connectivity, including delays. We also show that, in the suprathreshold regime, these networks are robust to heterogeneity in frequency and, in the subthreshold regime, they are robust to high interneuronal firing rates that occur if the interneurons participate in a large fraction of cycles. Next, we show that the cycle-skipping mechanism is built into the essential dynamics underlying type 2 excitability regardless of whether the neuron is biased in a subthreshold or suprathreshold regime. Finally, PV^{+} cells were much more likely to display both PIR and autapse-induced firing than GAD2^{+} cells, supporting the view that PV^{+} fast-firing basket cells are more likely to exhibit class 2 excitability than other types of inhibitory interneurons.

## Discussion

### Novelty

PIR can enable antiphase synchrony between two neurons (Perkel and Mulloney, 1974), as well as propagating activity patterns (Rinzel, Terman, Wang and Ermentrout, 1998). Wang and Rinzel (1993) showed that all-to-all networks of neurons with PIR due to a low-threshold calcium channel could synchronize via mutual inhibition. Heterogeneity resulted in clusters of silent versus active cells, with a slow variability in cluster composition. In contrast, firing patterns in our networks appear much more random and our networks do not require heterogeneity for sparse synchrony if the coupling is strong enough (Fig. 7*A*) because the quenched variability (van Vreeswijk and Sompolinsky, 1998) introduced by each neuron having a different set of presynaptic neurons suffices. Missed and skipped cycles (Glass and Winfree, 1984; Kaplan et al., 1996) are characteristic of the subcritical Hopf bifurcation (Ermentrout, 1996; Rinzel and Ermentrout, 1998; Prescott et al., 2008) underlying type 2 excitability (Hodgkin, 1948; Izhikevich, 2007). Fitzhugh (1976) described a “limiting trajectory” in a model with type 2 excitability, a concept that is very similar to the constant rebound interval after a strong inhibition observed in our study. The novelty of our theoretical work is to show that the tendency for missed cycles combined with the concept of a limiting trajectory enables robustness of sparse synchrony in an inhibitory network.

A damped oscillation (as in Fig. 1) in response to a small hyperpolarization is predicted for neurons near a subcritical Hopf bifurcation (Izhikevich, 2007; Rotstein, 2015), but we only found examples in the literature of damped oscillations (Gutfreund et al., 1995; Izhikevich et al., 2003) in neurons that exhibit spontaneous subthreshold oscillations, which do not characterize FS interneurons. PIR was first documented in hippocampal inhibitory interneurons by Cobb et al. (1995). The novelty of our experimental work is to demonstrate that the ability of inhibitory interneurons to exhibit both PIR and the damped oscillations characteristic of resonators (Fig. 1) allows participation in a network oscillation driven by PIR (Fig. 2), as well as a “limiting trajectory”—that is, constant rebound—in response to strong inhibition (Fig. 3).

### Historical perspective

Our interest in mechanisms of inhibitory interneuronal network gamma was sparked by a conundrum posed in a review of this literature by Bartos et al. (2007). They noted that neither the stochastic population oscillator (Brunel and Wang, 2003) nor the coupled oscillator (Wang and Buzsáki, 1996) model provided a satisfactory explanation of how inhibitory interneurons can participate on many cycles and remained synchronized to create a population rhythm despite documented variability in the level of external drive. The problem with a pure ING stochastic population model is that the SPC fired by a representative neuron is, on average, negatively correlated with the tightness of the window in which firing is disinhibited, whereas the tightness of synchrony is positively correlated with the same quantity. The negative correlation only holds true under the assumption that the firing times in each neuron approximate a Poisson process and inhibition in a network merely subtracts spikes during the periodic wave of inhibition. In the coupled oscillator model, each neuron is a rhythmic pacemaker. Synchronization due to interactions between neurons can then arise as an emergent property of the network and participation is no longer inversely correlated with the level of synchrony. However, the FS interneurons associated with *in vitro* models of ING have a steep frequency–current relationship, particularly at low levels of external drive and the physiological variability in drive is likely to be quite high. For example, in the context of the ING variant induced by bath application of metabotropic glutamate (mGLU), there is a 35% coefficient of variability (van Hooft et al., 2000) in the level of mGLU responses. Inhibitory networks of neurons with type 2 excitability successfully overcome both stated problems and may play a role in increasing the robustness of synchrony to high levels of both participation on a given cycle and heterogeneity between neurons.

### What current is responsible for the PIR-based mechanism?

Although currents like the hyperpolarization-activated nonspecific cation current (*I*_{h}) and the low threshold T-type Ca^{2+} current mediate PIR in some neurons, only postinhibitory rebound that is associated with currents responsible for spiking, as opposed to bursts or other slow oscillations (Wang and Rinzel, 1993; Rotstein and Nadim, 2014), is relevant for the robust cycle skipping mechanism proposed here. In the original Hodgkin–Huxley model, the relevant process is removal of inactivation of the fast, voltage-dependent Na^{+} channel. The Izhikevich model, unlike the classic Hodgkin–Huxley model, is not conductance based. However, the variable *u*, which functions as a net outward current, represents all the slower, restorative processes that oppose spike initiation, which include putative candidates for the PIR mechanism, including removal of inactivation of the sodium current.

### Network contributions from type 1 interneurons

There is a great diversity of inhibitory interneuronal subtypes (Klausberger and Somogyi, 2008) and there is likely diversity even within the PV^{+} FS basket cell population, especially because the excitability type can be modulated (Prescott et al., 2008). A recent study (Ferguson et al., 2013) intentionally modeled PV^{+} fast-firing C*A1* interneurons as integrators rather than resonators because they observed no evidence of subthreshold oscillations, unlike Figure 1, *A1* and *B1*. We also did not observe damped subthreshold oscillations unless the cell was depolarized from rest using a constant current before applying the virtual inhibitory conductance. An interesting direction for future research is to determine the effect on robustness if only a fraction of the inhibitory interneurons display type 2 excitability.

### Network contributions from pyramidal cells

In intact brain, excitatory pyramidal interneurons also participate in gamma oscillations. Introducing pyramidal neurons into stochastic population oscillator models (Brunel and Wang, 2003) allows the average level of excitatory drive to the interneurons to vary in time, which mitigates the conflict between participation and synchronization and can allow cycle skipping, but is limited by the degree to which excitatory neurons are allowed to synchronize their activity (Atallah and Scanziani, 2009).

One recent PING model (Economo and White, 2012) was constructed with pyramidal cells coupled to interneurons with type 2 excitability. Weakly synchronized pyramidal cells participated in ∼1 in 10 cycles, whereas interneurons demonstrated cycle skipping by participating with tight synchrony, on average, every other cycle, which is consistent with *in vitro* models of PING (Hájos et al., 2004). However, inhibition onto the interneurons was shunting and not hyperpolarizing. The present study focused on hyperpolarizing inhibition, so the mechanisms underlying cycle skipping in this network may be different, but may also arise from the intrinsic dynamics of interneurons with type 2 excitability.

The differential robustness of PING networks containing inhibitory interneurons with type 1 versus type 2 excitability has only begun to be explored. One study (Moca et al., 2014) using all-to-all coupled suprathreshold resonators embedded in PING circuits found improvement in the maintenance of a constant period in the presence of heterogeneity in extrinsic drive compared with networks with integrator interneurons. A second study of PING networks with heterogeneous drive to the pyramidal cells (Börgers and Walker, 2013) also found a more stable oscillation frequency and a more abrupt transition to suppression of the pyramidal cells as the external drive to the inhibitory interneurons was increased for resonator versus integrator interneurons. A third recent PING study (Baroni et al., 2014), like ours, emphasized the critical role of postinhibitory rebound and inhibitory, but not shunting, inhibition in increasing participation and synchrony in networks with resonator interneurons compared with integrators.

### Type 1 versus type 2 networks

There are two different ways (White et al., 1998) in which global synchrony is lost in an inhibitory network of type 1 oscillatory neurons with mild heterogeneity in frequency: phase dispersion and suppression of the slower neurons in the population by the faster ones. The robust cycle-skipping phenomenon described herein (Fig. 4) is not observed. Accompanying theoretical analyses (Chow et al., 1998) explained these results using a self-inhibited neuron as a reduced model of synchrony in inhibitory networks; relaxation back to rest was assumed to be exponential. These analyses need to be extended to include the damped oscillatory dynamics characteristic of resonators (Fig. 1). Moreover, the prevailing wisdom is that the population activity in networks of coupled oscillators does not change qualitatively in the presence of heterogeneity (Brunel and Wang, 2003). We show that the network oscillation can indeed change qualitatively, given type 2 dynamics and moderate noise. Specifically, the firing of individual neurons remains tightly locked to the population oscillation, but the neurons, on average, fire on fewer cycles as heterogeneity is increased. Moreover, we have shown the results are qualitatively similar in the subthreshold and suprathreshold regimes, also known as mean versus fluctuation-driven regimes.

### Model prediction

One important way to distinguish a stochastic population oscillation from the resonator network oscillation described here is suggested by Figure 10. We predict that single unit firing rate should vary inversely with vector strength of locking to the population at a constant frequency for a stochastic population oscillator, whereas a positive correlation between these quantities for individual type 2 neurons within a resonator network should characterize the oscillation mechanism presented herein.

### Conclusions

We present an elegant strategy for maintaining synchrony in the face of the myriad forms of biological inhomogeneity. The selective skipping of cycles in which the firing of a given neuron is not well aligned with the population allows the level of participation to adjust to the level of heterogeneity and noise. The most prominent drawback is the requirement for hyperpolarizing rather than shunting inhibition. This drawback may be mitigated if the membrane potential of PV^{+} interneurons is more depolarized *in vivo* than *in vitro*, given that a consistently small, but hyperpolarizing, driving force for ionotropic GABAergic synapses has been observed for interneurons in CA3 stratum lucidum *in vitro* across development (Banke and McBain, 2006). Moreover, the well documented existence of type 2 inhibitory interneurons suggests that the remarkable robustness of inhibitory networks of resonators may contribute to some of the many variants of gamma oscillations.

## Footnotes

↵*R.A.T.-H. and J.J.M. are co-first authors.

↵†J.A.W. and C.C.C. are co-last authors.

This work was supported by the National Institutes of Health (Grant R01NS054281 to C.C.C., Grant R01MH85074 to J.A.W., Grant R01MH085387 to C.C.C. and J.A.W., and Grant P30GM103340 to computational core).

The authors declare no competing financial interests.

- Correspondence should be addressed to Carmen C. Canavier, Department of Cell Biology and Anatomy, Louisiana State University, Health Sciences Center, 1901 Perdido Street, New Orleans, LA 70112. ccanav{at}lsuhsc.edu