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 jth neuron in a population of 300 neurons, time is in milliseconds, membrane potential vj is in millivolts, the applied current Ij is given in nanoamperes and is set to 0.15 unless otherwise noted, and the recovery variable uj is dimensionless. We have added a scale factor kj to some simulations to take into account that, in addition to variability in external drive, represented by Ij, 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 1K of Izhikevich, 2004, and Figure 8K of Izhikevich, 2007. vj is reset when the potential reaches threshold and the recovery variable uj 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 Ij 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 Jj(t) is a random process with zero mean and variance σN, the parameters for the synaptic connectivity term are Esyn for the synaptic reversal potential, gij for the synaptic conductance, sij(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. Esyn is set to −70 mV to simulate inhibition and the value of gij 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, gij = 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. gij 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/cm2, V is the cell membrane voltage in millivolts and t is time in milliseconds. The leak current is given by IL = gL(V − EL). The sodium current is given by INa = gNam∞3h(V − ENa). 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 IK = gKn4(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 ENa, EK, and EL were set to 55, −90, and −65 mV, respectively. The maximal sodium (gNa), potassium (gK), and leak (gL) conductances were set to 35, 9, and 0.1 mS/cm2, respectively. Ii 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 Esyn 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 8C, 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 Ij was set to give a 32–36 Hz spontaneous firing rate in isolated neurons of both model types (i.e., 5 μA/cm2 for integrator and 0.2 nA for resonator). Then the synaptic conductance gsyn was adjusted so that both networks had a frequency of 25–27 Hz (i.e., 1 nS/cm2 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 gsyn(vrest − Esyn) = 0.125 Ij for both models (i.e., gsyn = 0.125 μS/cm2 for integrator and 0.005 nS for resonator). In both cases, we then varied Ij 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 kj 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 10A, Ij was set below threshold to 0.15 nA for resonator and 1.5 μA/cm2 for integrator. Noise amplitude was set to give a mean firing rate of 10 Hz in both isolated model neurons (i.e., 3 μA/cm2 for integrator and 0.87 nA for resonator). Then, gsyn was set to produce the same level of synchronization for both networks (0.07 μS/cm2 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 gsyn for all synapses onto a given neuron from a Gaussian distribution truncated to eliminate negative values (Fig. 10A). For Figure 10B, Ij and gsyn had the same values as in Figure 10A, with gsyn set everywhere to its mean value, but the noise amplitude was swept through the range [2.8, 90] μA/cm2 for integrator and [0.06, 13] nA for resonator. We have summarized noise and connectivity parameters for all figures in Table 1.
Connectivity and noise parameters for all simulations
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 R2 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 R2 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 R2 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 NaH2PO4, 10 MgCl2, 25 NaHCO3, 12.5 glucose, 0.5 CaCl2), 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 NaH2PO4, 10 MgCl2, 25 NaHCO3, 25 glucose, 2 CaCl2) 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% O2/CO2 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 MgCl2, 0.2 EGTA, 10 HEPES, 20 KCl, 7 di(tris) phosphocreatine, 4 Na2ATP, 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 gmax is the maximal conductance, V is membrane voltage, Esyn 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 GABAA 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. 1A1). Moreover, an IPSP of sufficient amplitude produces a rebound spike (Fig. 1A2).
Resonant properties in response to IPSP. Dampened oscillatory response to a weak IPSP (gsyn = 0.1 nS) in the model (A1) and in an EC GAD2+ interneuron (B1). A stronger IPSP elicits a PIR spike in the model with gsyn = 0.5 nS (A2) and in an EC GABAergic interneuron with gsyn = 1.2 nS (B2). Inset shows feedforward circuit.
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 GABAA 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. 1B1) in many neurons. When the amplitude of the simulated conductance was increased, a rebound action potential was triggered (Fig. 1B2) 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. 2A1). 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. 2A1). Current noise was added to the model neuron to give a spread in the ISI histogram (Fig. 2B1). We then repeated this protocol in a biological setting, again using the virtual GABAA 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. 2A2). Indeed, the histogram of ISIs observed in the self-connected neuron was narrow (Fig. 2B2), 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.
Simulated network feedback using a virtual, spike-triggered IPSP. A, Membrane potential (top trace) and injected conductance waveform in the model (σN = 0.1 nA, gsyn = 0.5 nS) (A1) and in a representative EC GAD2+ interneuron, gsyn = 5 nS (A2). Inset shows simulated virtual autapse. Triangles indicate an initial input not triggered by a spike. B, ISI probability histogram for the model (B1) and a different EC GAD2+ (B2) neuron from the one shown in the previous figure.
Using a model neuron with an autapse as in Figure 2A, an additional, untriggered input was applied after a feedback IPSP, but before the corresponding rebound spike was emitted (triangle in Fig. 3A). 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. 3A2). 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.
External input resets the of phase self-sustained oscillations. A, Model. A1, Each spike triggers a feedback IPSP (gsyn = 0.5 nS) as in Figure 2. An additional, untriggered synaptic input (triangle, gsyn = 1.0 nS) resets the phase of the autapse-induced oscillation in a model neuron (σN = 0.16 nA). A2, The phase-resetting curve measured for a noisy model neuron with an autapse is linear with a slope near 1, corresponding to a constant latency after the untriggered input. The best fit line is in blue and the potential waveform is in red. B, Experiment with similar parameters as the model. Each spike triggers a feedback IPSP (gsyn = 0.6 nS), as in Figure 3A1, and an additional, untriggered synaptic input (gsyn = 1.2 nS) resets the phase of the autapse-induced oscillation. B1, Measurement protocol for the PRC in an EC PV+ interneuron. B2, PRC for the PV+ neuron is also linear. One outlying point in B2 was omitted.
An additional, untriggered input was applied (triangle in Fig. 3B) 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 gsyn = 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.
Robust synchronization in noisy sparsely connected network due to spike skipping. All neurons fire a PIR spike on every cycle in the absence of noise. When noise is introduced, some random subset of population fires after each population spike, which is evidenced in the single neuron trace (A1), the network raster plot (A2), and spike rate histogram (A3). B, The histogram of all ISIs in all neurons within the population shows a clear peak at the network frequency. Peaks at multiples of this frequency indicate that one or more cycles were skipped.
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. 5A1). The trace of the membrane potential in a representative neuron (Fig. 4A1), as well as the raster plot (Fig. 4A2), show that cycles are skipped. The spike time histogram (Fig. 4A3) 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. 4B) 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.
Inhibition selectively removes spikes that occur later in the cycle. A, No delay. B, Conduction delay of 3 ms. SD of noise = 0.375 nA. Connectivity is sparse with 40 connections per interneuron.
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. 5A). 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. 5B) 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. 6A), 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 v2 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.
Cycle-skipping mechanism. A, Subthreshold. Trajectories leading to and from the saddle (filled diamond) determine whether trajectories of individual neurons fire or skip a cycle. B, Bistable. The unstable limit cycle (partially dashed/partially solid closed curve) intersects with the voltage nullcline (solid gray) where the net current = v′ = 0. The dark part of both curves forms the boundary for cycle skipping. The dashed line shows u = bV, where u is at its steady state u′ = 0. The intersections of v′ = 0 and u′ = 0 give the fixed points. C, Suprathreshold. Again, trajectories leading to and from the saddle (filled diamond) determine whether trajectories of individual neurons fire or skip a cycle. a.u., Arbitrary units.
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. 6B); 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. 6C), 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. 7A). 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 gsyn 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 R2 > 0.9 and a region of intermittent synchrony in which R2 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.
Effect of synaptic parameters on synchronization in a sparsely connected noiseless network. A, Parameter space divides into regions of sparse synchrony and intermittent sparse synchrony and silence. B, ISI histograms. B1, Global synchrony. B2, Cycle skipping. B3, Cycle skipping with a broader window. B4, Broad peaks superimposed on a baseline level of firing reflect intermittent sparse synchrony.
For gsyn, there a minimum value below which PIR is not generated (Fig. 1) but it is barely visible on the scale in this figure. As gsyn 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. 7B). Peaks other than the first one indicate skipping one or more cycles and broader peaks reflect a broadening of the window. Increasing gsyn increases the degree of quenched randomness and destroys the global attractiveness of full synchrony observed in the noiseless part of Figure 4A. For an instantaneous conductance (delta function), there is no upper bound on gsyn because the effect of the synaptic conductance saturates as the membrane potential approaches the reversal potential. However, synchrony is weakened as gsyn 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 7A, instead sparse synchrony waxes and wanes, with less synchrony as gsyn continues to increase. This effect of increasing gsyn 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 6A 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 4A. 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. 8A). The R2 vector strength (solid black curve) measure of synchronization remained above 0.7 even for this strong level of noise. Next (Fig. 8B), 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 8C, 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 R2 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 1B). Nonetheless, the SPC decreased as heterogeneity in connectivity was increased, but the synchronization remained near perfect for a CV of 1. Finally (Fig. 8D), 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.
Response of the network oscillation to heterogeneity. Only one type of heterogeneity was introduced in any given simulation. Vector strength R2 is a measure of synchronization shown by the black curves and average spikes of individual neurons per cycle (SPC) of the population oscillation (dashed curve). A, Noise. B, External drive. C, Synaptic conductance. D, Conduction delays. In B, the R2 is variable at low SPC and depends upon initial conditions.
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 8B, there was a slight increase in R2 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 R2 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 R2 plummeted but participation (SPC) remains high, and suppression, in which R2 remains high but the faster neurons suppress other neurons so the participation drops. The R2 vector strength for a network of integrators plummeted near 0 with a CV of only ∼0.1 or 10% (dashed curve in Fig. 9A). Conversely, a comparable tuned network of resonators (see Materials and Methods) was able to maintain strong synchronization (solid line near 1 in Fig. 9A) at a CV of 0.8, at a cost of decreasing the SPC (Fig. 9B). 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. 9B and C2 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.
Robustness of synchrony in coupled oscillator networks of resonators versus integrators (all-to-all connectivity). Increasing the CV of the intrinsic frequencies affected synchronization (A) and participation (B) differently. Raster plots (C) show how synchrony fails for integrators in the tonic regime (C1), integrators in the phasic regime (C2), and resonators (C3). Note that only one set of curves is plotted for the resonator networks, although two parameter regimes were explored as described in the Materials and Methods because the curves were indistinguishable.
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 R2 vector strength of 0.97. At this value of R2, 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 10A) compared with the resonator network (solid line in Fig. 10A) 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 R2. 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 R2 vector strength (dashed curve in Fig. 10B). 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 R2, but only slightly, until an inflection point was reached. After the inflection point, R2 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 kj 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 10B decreased approximately proportionately (data not shown).
Robustness of synchrony in stochastic network oscillator of resonators versus integrators. These networks have all-to-all rather than sparse connectivity. In a subthreshold regime, the CV of the conductance (A) and the SPC (B) were systematically varied. C, Peaks in the ISI histogram became broader and the valleys filled in as noise was increased.
It is clear in the histograms in Figures 7B and 10C 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 6C, 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. 7A) 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 (Ih) and the low threshold T-type Ca2+ 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 CA1 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
