Mechanisms of Synchrony in Heterogeneous Inhibitory Interneuronal Networks for Type 1 versus Type 2 Excitability

PV+ fast spiking basket interneurons are often implicated in gamma rhythms. Here we focus on mechanisms present in purely inhibitory networks. Neurons with type 1 excitability can fire arbitrarily slowly, whereas those with type 2 excitability cannot fire below a minimum frequency. We systematically examine how excitability type affects synchronization of individual spikes to a population rhythm in the presence of heterogeneity and noise, using model neurons of each type with matched F/I curve, input resistance, time constant and action potential shape. Population synchrony in noisy heterogeneous networks is maintained because neurons either fire within a tight time window or skip that cycle. Type 2 neurons with hyperpolarizing inhibition skip cycles due to their intrinsic dynamics; we show here the cycle skipping mechanism for type 1 neurons or type 2 neurons with shunting inhibition is synaptic and not intrinsic. Type 2 neurons are more resistant than type 1 to partial and complete suppression in networks with hyperpolarizing inhibition that exhibit network gamma. Moreover, type 2 neurons are recruited more rapidly and more completely into theta-nested gamma. In contrast, type 1 networks perform better with shunting inhibition on both counts, because the nonlinear dynamics in that case favor suppression of type 2 compared to type 1 neurons. Conductances that control excitability type may provide a therapeutic target to improve spatial and working memory and other tasks that rely on gamma synchrony or phase amplitude coupling. Author Summary The collective, synchronized activity of neurons produces brain rhythms. These rhythms are thought to subserve cognitive functions such as attention and memory encoding and retrieval. Faster rhythms are nested in slower rhythms as a putative way of chunking information. A subset of neurons called fast spiking basket cells tend to inhibit other neurons from firing. These neurons play an important role in oscillations, and in the coupling of faster oscillations to slower ones. In some brain regions these neurons can fire arbitrarily slowly (type 1 dynamics) whereas in others they cannot fire below a minimum cutoff frequency (type 2 dynamics). Mathematically, these distinct origins of rhythmic firing are signatures of very different dynamics. Here, we show that these distinct excitability types affect the ability of networks of these neurons to synchronize their fast oscillatory activity, as well as the ability of slower oscillations to modulate these fast oscillations. The exact nature of the inhibitory coupling, which may vary between brain regions, determines which type synchronizes better and is modulated better.

Abstract PV+ fast spiking basket interneurons are often implicated in gamma rhythms. Here we focus on mechanisms present in purely inhibitory networks. Neurons with type 1 excitability can fire arbitrarily slowly, whereas those with type 2 excitability cannot fire below a minimum frequency. We systematically examine how excitability type affects synchronization of individual spikes to a population rhythm in the presence of heterogeneity and noise, using model neurons of each type with matched F/I curve, input resistance, time constant and action potential shape. Population synchrony in noisy heterogeneous networks is maintained because neurons either fire within a tight time window or skip that cycle. Type 2 neurons with hyperpolarizing inhibition skip cycles due to their intrinsic dynamics; we show here the cycle skipping mechanism for type 1 neurons or type 2 neurons with shunting inhibition is synaptic and not intrinsic. Type 2 neurons are more resistant than type 1 to partial and complete suppression in networks with hyperpolarizing inhibition that exhibit network gamma. Moreover, type 2 neurons are recruited more rapidly and more completely into theta-nested gamma. In contrast, type 1 networks perform better with shunting inhibition on both counts, because the nonlinear dynamics in that case favor suppression of type 2 compared to type 1 neurons. Conductances that control excitability type may provide a therapeutic target to improve spatial and working memory and other tasks that rely on gamma synchrony or phase amplitude coupling.

Author Summary
The collective, synchronized activity of neurons produces brain rhythms. These rhythms are thought to subserve cognitive functions such as attention and memory encoding and retrieval. Faster rhythms are nested in slower rhythms as a putative way of chunking information. A subset of neurons called fast spiking basket cells tend to inhibit other neurons from firing. These neurons play an important role in oscillations, and in the coupling of faster oscillations to slower ones. In some brain regions these neurons can fire arbitrarily slowly (type 1 dynamics) whereas in others they cannot fire below a minimum cutoff frequency (type 2 dynamics). Mathematically, these distinct origins of rhythmic firing are signatures of very different dynamics. Here, we show that these distinct excitability types affect the ability of networks of these neurons to synchronize their fast oscillatory activity, as well as the ability of slower oscillations to modulate these fast oscillations. The exact nature of the inhibitory coupling, which may vary between brain regions, determines which type synchronizes better and is modulated better.
Gamma frequency  oscillations are thought to serve as substrates for working memory, conceptual categorization, and attention [9], and are altered in psychiatric disorders, for example schizophrenia, dementia, and autism [10]. PV+ FS basket cells play an important role in theta nested gamma [11]. Moreover, nesting of gamma within theta has been proposed a substrate for episodic memory [12] and disruption of theta nested gamma has been proposed to explain deficits in spatial memory in temporal lobe epilepsy [13,14].
Here we focus on mechanisms present in purely inhibitory networks, and the role of excitability type of the FS interneurons in gamma frequency synchronization and cross frequency synchronization with theta oscillations. Neurons with Hodgkin's type 1 excitability are able to spike arbitrarily slowly, whereas those with type 2 excitability have an abrupt onset of repetitive firing that cannot be maintained below a threshold frequency [15,16]. Our recent work [17] shows that FS neurons likely exhibit type 2 excitability in the medial entorhinal cortex, and others have shown that they exhibit type 2 excitability in striatum [18] and neocortex [19]. However, previous work supports type 1 excitability in hippocampal area CA1 [20][21][22]. We will focus on reduced, 2D models here.
The advantage of 2D models is that they can be visualized in a two-dimensional phase plane that reveals the essential features of the intrinsic dynamics. For type 2 excitability, there is a circular flow near the associated subcritical Hopf bifurcation that results in subthreshold oscillations, postinhibitory rebound inherent in the spike generating mechanism, and a resonant spiking frequency.
These features are absent in neurons with type 1 excitability. In contrast, type 1 neurons rely on a saddle-node-on-an-invariant-circle (SNIC) bifurcation, which creates a slow channel in the phase plane where the net current is very small, allowing the neurons to spike arbitrarily slowly.
Previously [17], we noted that cycle skipping was a way for a network to preserve synchrony of individual spikes with the population in the presence of heterogeneity and noise. In other words, spikes that would have occurred too late to be in synchrony with the population are suppressed.
In the case of type 2 excitability, this can be accomplished via the intrinsic dynamics, because in all regimes (excitable, bistable and oscillatory) there is a boundary between trajectories that lead to spiking and those that lead to cycle skipping. However, we did not calibrate the models we compared in such a way as to maximize the fairness of the comparison. Moreover, we only looked at hyperpolarizing and not shunting inhibition.
Fast, ionotropic inhibition in the central nervous system in generally mediated by GABAA receptors, with chloride ions as the charge carrier. The reversal potential of these channels depends on the intracellular concentration of chloride. For quiescent neurons, for hyperpolarizing inhibition the synaptic reversal potential that is negative to the resting potential, whereas for shunting inhibition the reversal potential is more positive than the resting potential but negative to the spike threshold. For oscillatory neurons, the membrane potential during the interspike interval substitutes for the resting potential [23]. Shunting inhibition is sometimes defined as an increase in synaptic conductance in the absence of an obvious change in membrane potential [24]. In our model neurons, a reversal potential of -75 mV produces hyperpolarization whereas -65 mV does not, so we classify the latter as shunting in our context. The reversal potential for chloride in fast spiking basket cells in vitro is about -52 mV [25] in the dentate gyrus, which implies shunting inhibition in those cells. Previous studies in CA1 and CA3 concluded that inhibition between basket cells was hyperpolarizing [26,27]. However, the Cl-concentration can be modulated [28], and there can be intracellular gradients in chloride concentration [29]. Therefore, there is sufficient uncertainty regarding the precise reversal potential of chloride at synapses between interneurons to warrant a systematic study of both types of inhibition.
Here, we systematically examine how excitability type affects synchronization of individual spikes to a population rhythm in the presence of heterogeneity and noise. Simple 2D model neurons of each type were calibrated to have a very similar F/I curve, input resistance, time constant and action potential shape in order to isolate the consequences of excitability type alone on the robustness of synchronization in inhibitory interneuronal networks.

Phenomenological model of type 1 and type 2 excitability
We began with the standard 2D reduction [30] of the 4D Hodgkin-Huxley model [31]. The activation of the voltage-dependent sodium current was assumed to be fast and set to its steady-state value with respect to the membrane potential ( . The inactivation variable for the voltagedependent Na + current ( h ) was yoked to the variable ( n ) for the activation of the delayed rectifier K + current, under the assumption that the slow time scale for these variables was similar. This 2D system with one fast variable (membrane potential, v ) and one slow variable ( n ) is amenable to phase plane analysis under fast/slow assumptions [32]: Figure 1 shows a phase plane analysis of the bifurcations underlying type 1 (Fig. 1A) and type 2 ( Fig 1B) excitability. In the phase plane, fixed points occur at the intersection of the n-nullcline (gray curve), which is simply the steady state activation curve for the n activation variable, and the voltage nullcline at rest (red curve), which is the set of values of n and v for which the net ionic current plus any applied current is zero. The voltage nullcline is N-shaped with three branches, but in order to emphasize the bifurcation point, only the stable left and unstable middle branch and their associated fixed points are shown in Fig 1. For this fast/slow system the leftmost branch is stable, and the middle branch is unstable due to the autocatalytic sodium current. At rest (with no applied current), the fixed point (filled circle) on the left branch is stable and determines the resting potential at -68 mV for both the type 1 (Fig. 1A1) and type 2 ( Fig. 1B1) cases.
The first important difference is that for type 2 there is a single fixed point at all values of applied current, whereas for type 1 the number of fixed points depends upon the level of the applied current.
For type 1 at rest (Fig. 1A1), in addition to the stable fixed point that determines the resting potential, there is an unstable fixed (open circle) point on the middle branch, as well as another one on the unstable branch that is not shown. As the applied current is increased to the stable and unstable fixed points shown collide (Fig. 1A2) and form a saddle node (circle half filled). This collision coincides with the birth of limit cycle surrounding the rightmost unstable fixed point on the unstable branch (not shown), therefore the onset of spiking occurs at a saddle-nodeon-an-invariant-limit -cycle (SNIC) bifurcation [16,33]. As the applied current is increased further and the colliding fixed points annihilate each other. The neuron is left with no stable resting potential. The intersection of the limit cycle with the saddle node gives rise to a trajectory with infinite period. As the applied current is increased further, the trajectory must pass through a gap between the two nullclines; the rate of passage though this bottleneck is arbitrarily slow with increasing proximity to the bifurcation, hence the arbitrarily slow frequencies obtainable by type 1 neurons shown in Fig. 2A (black dots).
In contrast, for type 2 excitability the Hopf bifurcation occurs when the lone stable fixed point corresponding to the rest potential loses stability (Fig. 1B1). When the applied depolarizing current , the fixed point (Fig. 1B2, open circle) moves past the trough of the nullcline onto the unstable branch of the voltage nullcline. The Hopf bifurcation is subcritical [16,33,34] because in the range of input currents between   2 1.74, 2.11 / I A cm   the model exhibits bistability between stable limit cycle corresponding to repetitive spiking and a stable fixed point ( Fig. 2A, overlap in the red cross at 0 frequency and at a nonzero frequency). These attractors are separated by unstable limit cycle (not shown). The stable and unstable limit cycles collide and annihilate each other in a saddle node of periodics at . Excitability type is defined as the response of a quiescent neuron as the applied current is increased [15]; the sudden onset of spiking as the applied current is increased occurs at about 30 Hz when the quiescent state loses stability as shown in Fig. 1B.

Steady State Synchrony in Oscillatory Networks
We constructed networks of 300 sparsely and randomly connected neurons with heterogeneity in frequency by distributing the bias current uniformly along a region of the F/I curve that spanned a 20 Hz range (see black bars in Fig. 2A). The conduction delays between neurons were also uniformly distributed between 0.7 and 3.5 ms. In our previous work [35] in a homogeneous network, delays on the short end of this range favor a solution with two subclusters in antiphase, whereas delays at the longer end of the range favor global synchrony of a single cluster. Using a mixture of delays results in solutions that are not obviously one or two clusters, but are transitional between these two extremes.  Fig. 3A and B) in is much flatter for type 1 than type 2, because the type 2 histogram is skewed to higher participation rates.

Steady State Synchrony in Oscillatory Networks-Hyperpolarizing Inhibition
This tendency for greater suppression of type 1 neurons with hyperpolarizing inhibition was pre-  Fig. 1A1) at the intersection with the n-nullcline. The unstable branch of the V-nullcline (arrow indicating firing threshold points to this branch) separates neurons into two groups. Neurons whose trajectory has already moved to the right of this branch will fire an action potential, but those that fall on the left side of this boundary when the channel closes will skip this network cycle. In general, the neurons that skip are the slower neurons; they will move leftward towards the stable fixed point. The two groups are clearly shown in For example, there is no slow channel. Instead, as explained in Fig. 2A, near a subcritical Hopf bifurcation, a bistable region exists in which quiescence at a stable fixed point (arrow) co-exists with repetitive spiking at the same value of net applied current. Whether the neuron is silent or active depends on recent history, meaning the current location of the trajectory. Specifically, the red closed curve in Fig. 4B1 is an unstable limit cycle that divides the neurons into two groups.
Those inside the red curve will spiral into the stable fixed point at its center, whereas those outside will curve around it to fire an action potential. Thus the mechanism for cycle skipping has a component that results from the intrinsic dynamics. There is no bias toward faster cells, all neurons are almost equally likely to fall outside of the quiescent zone and fire an action potential. Fig. 4B2 shows that as neurons move to the right and fire actions potentials, inhibition accumulates as in

Steady State Synchrony in Oscillatory Networks-Shunting Inhibition
The results in the previous section were for hyperpolarizing inhibition. For both type 1 and type 2, the bifurcation that gives rise to spiking occurs at about -64 mV. Hyperpolarizing inhibition with a reversal potential of -75 mV hyperpolarizes the membrane at most points during the interspike interval except during the trough of the AHP. Shunting inhibition with a reversal potential of -65 mV does not produce big changes in the membrane potential.  4B2 is that the quasithreshold described above (red curve in Fig. 6B2) is now visible as the "ghost" of the unstable limit cycle at the inhibition midpoint. Moreover, the inset in Fig. 6C2  neurons the cutoff frequency is also variable. We incorporated this additional dimension of variability into our simulations using a scale factor for the model dynamics as described in the Methods.
In order to determine the effect of heterogeneity in this parameter, we scaled the dynamics of both variables for each neuron in order to achieve uniformly distributed cutoff frequencies in the range from 10 to 60Hz for type 2 neurons, and used the same range of scale factors for type 1 neurons. We then redid the heatmaps for each type with hyperpolarizing and with shunting inhibition (supplementary figures 1 and 2). Although both types show a reduction in synchronization, the qualitative results that type 2 networks are more resistant to partial and complete suppression for hyperpolarizing but not shunting inhibition still hold.
We also checked to make sure that a saddle-node-not-on-an-invariant-circle (SN) would have qualitatively similar synchronization tendencies. We adjusted the parameters ( 0 3. ms ) of the single type 1 model neuron such that the stable spiking limit cycle is not tangent to the saddle node at the bifurcation. Nonetheless, at synaptic conductance strengths sufficient to create a stable fixed point during network activity, a "bottleneck" or slow channel emerges with the same functional role in organizing network activity as if there were a saddle node on an invariant circle (not shown).

Phase Amplitude Coupling
In the hippocampus, gamma power is maximal when nested in theta oscillations [37]. In order to determine the relative abilities of networks of neurons with type 1 excitability versus type 2 excitability to produce theta-nested gamma, we drove these networks with perfectly sinusoidal inhibitory waveforms at a fixed frequency in the theta range, in the presence of the constant heterogeneous depolarizing bias currents distributed as described in the Methods. The depolarization mimics tonic activation of metabotropic glutamatergic/cholinergic receptors. PV+ basket cells in freely moving rats fire at about 7 Hz during low oscillatory periods, but that rate triples to 21 Hz during theta oscillations [38], presumably due to greater tonic excitation. The sinusoidal drive mimics phasic inhibition from the septum. The amplitude of modulation is set to be the same for each neuron within the population, to reflect the common source of the modulation. The noise and conductance parameters were selected in a regime in which all four types of networks synchronized well. The attributes of the network oscillations for the selected parameter regime are given in Table 3. In additional sets of simulations (not shown), we confirmed that the results shown below are qualitatively similar for wide range of synaptic conductance and noise levels and are not specific for the chosen values. Figure 7C shows the results in the two dimensional parameter space of modulation depth and modulation frequency, with the color in panel C1

Phase Amplitude Coupling-Hyperpolarizing Inhibition
and C2 indicating the un-normalized vector strength as described in the methods. If the vectors were normalized to reflect only how tightly locked the LFP envelope was to the theta drive, the two types perform equally. Removing the normalization reveals the greater recruitment of the population into the nested gamma in type 2 networks. The heatmaps in Fig. 7C1 and C2 show that low frequencies and shallow modulations are most effective, especially for type 2. The heatmap in Fig. 7D show that Type 2 always outperforms type 1 because the red color is always greater than zero.
Almost no theta modulation is evident in the membrane potential traces of individual neurons Fig   7A and B. As the network recovers from inhibitory modulation, the firing rate in the population increases, therefore network inhibition fills in when the external inhibitory drive wanes. This feedback mechanism is also responsible for keeping neurons near the bifurcation at minimal inhibition in steady-state oscillations (see Figs. 4A1 and B1 and 6A1 and B1).
In order to explain the superior performance of type 2 for hyperpolarizing inhibition, we can refer back to the phase plane in Fig. 4A. In order to synchronize, type 1 relies on the accumulation of inhibition as neurons escape through the slow channel, whereas there is no slow channel for type 2 (Fig. 4B). Instead, all neurons have access back to the strong vector fields away from the fixed point and an opportunity to spike. We hypothesize that the improved recruitment of neurons into the gamma oscillation in type 2 networks was because type 1 networks need more time to establish synchrony. We tested this idea in a network in which all neurons had the same applied bias current and therefore the same intrinsic frequency (see Figure 8A). The speed of synchronization is clearer for the homogeneous frequency case because participation is more uniform. A headto-head comparison over two cycles of sinusoidal revel reveals that indeed synchronization begins earlier in the type 2 homogeneous network compared to the type 1 homogeneous network (see leftmost vertical gray bar), which accounts in part for the larger vector strength for theta phase gamma amplitude coupling for type 2 with hyperpolarizing inhibition. The second gray bar shows that the synchrony persists longer for type 2 as well.
The superiority of type 2 for hyperpolarizing inhibition derives from the ability to more quickly recruit the activity of sufficient neurons to establish synchrony at gamma frequencies. The presumption is that the interneurons are in an excited, oscillatory state during theta activity, and are rhythmically inhibited by the septum. For modulation that is too shallow, gamma activity is ongoing and theta power is weak (bottom blue strips on Fig.7 C1 and C2). For modulation that is sufficiently deep to establish theta power, the modulation must be sufficiently slow or shallow that there are still long enough windows of time above the spiking threshold to recruit enough active neurons to establish gamma synchrony (red areas in bottom left half of Fig. 7 C1 and C2).

Phase Amplitude Coupling-Shunting Inhibition
The faster recruitment of type 2 oscillators did not persist for homogeneous networks with shunting inhibition (Fig. 8B). Instead, type 1 oscillators are recruited more quickly into nested theta gamma because shunting inhibition clamps neurons closer to the spiking threshold (see Fig. 6C1 vs. Fig. 4C1). Type 2 networks clearly lose their superiority, likely because of the tendency of type 2 oscillators with shunting inhibition to remain trapped at the fixed point inside the unstable limit cycle ( Figure 6B1 and C2). Figure 9A and B show examples of Type 1 and 2 networks with shunting inhibition. Just as switching from hyperpolarizing to shunting inhibition nullified the advantage of type network in robustness to suppression during steady state synchrony shown in Figures 3 and 5, changing the inhibition from hyperpolarizing to shunting nullifies the advantage of type 2 for phase amplitude coupling as predicted by Fig. 8. In fact, Fig. 9D shows that type 1 networks have better phase amplitude coupling for deep modulations (blue indicates a type 2 -type 1 difference is less than zero). Although type 2 appears to perform better at low modulation amplitudes, the vector strength is so small (blue at the bottom of Fig. 9C1 and C2) as to render this regime unimportant. The tendency of the unstable limit cycle in Fig. 6B1 to trap trajectories clearly reduces the speed of synchronization of type 2, and is likely responsible for its poorer performance with shunting inhibition. Supplementary figures 3 and 4 shows that the PAC tendencies for type 1 versus type 2 are preserved under the second type of heterogeneity utilized in supplementary figures 1 and 2.

Summary of Results
We present two major results. The first is that networks of heterogeneous neural oscillators with type 2 excitability are more resistant to suppression of the slower oscillators than those with type 1 excitability when they are coupled with hyperpolarizing inhibition but not with shunting inhibition. The second result is that theta phase to gamma amplitude coupling is more strongly recruited in the type 2 networks, again for hyperpolarizing but not shunting inhibition. Moreover, we provide mechanistic explanations for these phenomena, which frame conditions for the generality of these results.
The tendency of slower type 1 neurons to be more suppressed than type 2 neurons for hyperpolarizing inhibition relies on two principles. The first is the creation by inhibition of a slow channel near the SNIC bifurcation for type 1. The channel serves to line the fixed points of the heterogeneous oscillators up along the slow nullcline and along the direction of motion of the fast subsystem. This arrangement clearly favors the faster oscillators, and is less prominent for shunting versus hyperpolarizing inhibition.
The second principle is the circular motion about the subcritical Hopf bifurcation for type 2 neurons and the existence of a boundary in the phase plane between spiking and skipping neurons allows for broader participation in steady state network oscillation for hyperpolarizing but not shunting inhibition. The critical assumption is that the stable fixed point just prior to the Hopf bifurcation lies near the shunt reversal potential, but well above the hyperpolarizing reversal potential. This assumption seems reasonable based on the definitions of hyperpolarizing and shunting inhibition combined with the location of the Hopf at the point that positive feedback from the sodium channel destabilizes the rest potential. Thus hyperpolarizing inhibition pushes all neurons leftward (in the fast direction) and the circular motion brings them downward and into the fast vector fields that give neurons a much more equal chance to fire than type 1. However, shunting inhibition, rather than freeing trajectories from the unstable limit cycle surround the stable fixed point just prior to the Hopf, traps them within that limit cycle, favoring suppression.
These principles apply to the transient synchronization in theta-nested gamma as follows. Note that the steady state synchronization and participation are similar between type 1 and type 2.
However, type 2 neurons are more quickly recruited into nested theta gamma for hyperpolarizing inhibition because their intrinsic dynamics do not require accumulation of synaptic inhibition to split the populations into spiking and skipping groups. Faster recruitment also allows synchronization, which requires a minimum number of active neurons, to persist longer. The phase plane portrait also accounts for the slower recruitment of type 2 neurons with shunting inhibition into nested theta gamma and the shorter persistence that accompanies late recruitment. The mechanism is again the tendency of shunting inhibition to clamp the type 2 trajectories inside the unstable limit cycle, which slows recruitment into theta nested gamma.

Generality
A heterogeneous, multicompartmental neuron can be considered as a chain of diffusively coupled oscillators. Although the SNIC bifurcation is generic in single compartment models, studies of heterogeneous oscillators, each exhibiting a SNIC in isolation [39], exhibit a complicated solution structure in a diffusively coupled system. It is not clear that a heterogeneous, multicompartmental neuron can exhibit a SNIC precisely, so it was important to show that the principles described in the preceding section for type 1 neurons still hold for a network in which individual neurons do not have a SNIC, but rather a SN not on an invariant limit cycle. The slow channel still forms.
Another important way in which we checked the generality of the results was to introduce a second type of heterogeneity. We varied the slope of the f/I curve for type 1 neurons and both the cutoff frequency and slope for the type2 neurons by rescaling the temporal dynamics. The cutoff frequency for fast spiking PV+ basket cells [40] had an estimated standard deviation of 40%.
Thus, the 10-60 Hz range in cutoff frequencies explored should map onto biologically plausible distributions of cutoff frequencies. Qualitatively, our major results were preserved under these manipulations ( Supplementary Figs. 1-4).

Previous comparisons of type 1 versus type 2 excitability
Previously Rinzel and Ermentrout [41] contrasted the phase portraits and bifurcation structure for type 1 versus type 2 excitability using the Morris Lecar [42] model in the two regimes. Izhikevich [16] also showed in the phase plane that various minimal conductance-based models could exhibit saddle node and Hopf bifurcations. Neither study explored the implications of the bifurcation type for synchronization.
In the weak coupling regime, Ermentrout and colleages [43,44] found that type 2 neurons receiving noisy common input synchronize better than type 1. Borgers and Walker [45] compared the synchronization tendencies of type 1 [20] and type 2 [46] interneurons embedded in excitatory/in-hibitory networks with type 1 excitatory neurons. They found that when gap junctions were incorporated into their model, there was a sharper transition from pyramidal interneuronal network gamma to interneuronal network gamma in the type 2 networks compared to type 1 networks.
Rich et al [47] found that excitatory/inhibitory networks with type 2 interneurons were more robust to changes in network connectivity compared to type 1. None of the type 1 versus type 2 model comparisons cited above were calibrated to match spike shape, shape of the F/I curve, time constant and input resistance across excitability types as were the models in our study. Thus, we extend these previous comparisons so that any difference in network activity is due to the bifurcation type.

Cycle skipping mechanisms
Cycle skipping to enforce population synchrony within the inhibitory neuron population has been observed in several studies in networks with both type 1 [48], and type 2 interneurons [17,49].
However, the bifurcation structure for the four cases (combinations of type 1 and type 2 excitability with either hyperpolarizing or shunting inhibition) has not been examined in detail previously. It is clear that interneurons do not fire on every gamma cycle in vivo. For example, CA3 FS cells fired at 21 Hz during 45 Hz gamma at 45 Hz, thus they fired roughly on every other cycle [50].
Cycle skipping provides a mechanism by which coupled oscillator models can produce tightly synchronized firing with sparse firing in individual neurons. This mechanism is distinct from a clustering mechanism that was advanced to introduce robustness of population synchrony. In that mechanism, interneurons did not participate on every gamma cycle because synchronized subclusters fired in sequence [51]. Previously [17], we found that cycle skipping more prominent in Type 2 networks, but it was not possible to do a fair comparision across the parameter space with the unmatched models we used in that study. The appropriate conclusion is that cycle skipping is intrinsic to the dynamics of type 2 neurons connected by hyperpolarizing inhibition, but in all other cases synaptic mechanisms predominate over intrinsic mechanisms. Cycle skipping allows robustness to heterogeneity in excitatory drive and connectivity for both type 1 and type 2 networks, overcoming a perceived weakness of coupled oscillator models of network gamma since they were introduced (see next section).

Previous studies on suppression
Computational studies have often shown that fast oscillations based on reciprocal inhibition are exquisitely sensitive to heterogeneities in the network [1,20,23,36]. Inhibitory neural networks generally lose synchrony as heterogeneity is increased in one of two ways, phase dispersion or suppression, depending on the ratio of the time constant for decay of inhibition to the population frequency [36,52] . Here, we looked at fast synapses, which tend to favor the suppression regime over phase dispersion. Early work [20,36,52] on robustness of inhibitory interneuronal synchrony focused of type 1 excitability only. These studies did not include conduction delays, and concluded that hyperpolarizing inhibition and non-physiological low levels of heterogeneity in excitatory drive were required for synchrony. Subsequent studies [25] used the same type 1 Wang and Buzsaki model added a ring topography with conduction delays ranged from 0.7 to 10.5 ms; they found that networks with shunting inhibition were more robust to heterogeneity than with hyperpolarizing inhibition. Their raster plots clearly show more suppression for type 1 for hyperpolarizing compared to shunting inhibition, in agreement with our findings. They repeated their simulation with the type 2 model of Erisir et. al. [46], and obtained similar synchronization results, as in our study.
However, they did not show raster plots or report on suppression in the type 2 networks. We also used a distribution of delays (from 0.7 to 3.5) to prevent the formation of two cluster solutions [35] and stabilize global synchrony by moving the operating point away from the destabilizing discontinuity in the phase resetting curve at 0 and 1.

Therapeutic directions for PAC and cognition
The hippocampal theta rhythm has been shown to be necessary for spatial learning by rats in a water maze [53]. Phase amplitude coupling between theta and gamma likely plays an important role in cognition [54], and we have shown here that type 2 excitability in interneurons strengthens phase amplitude coupling between theta and gamma rhythms for hyperpolarizing but weakens it for shunting inhibition. This suggests the possibility that interneurons could modulate their excitability type according to whether inhibition is hyperpolarizing or shunting, or vice versa. Moreover, this suggests a therapeutic strategy to manipulate excitability type, and thereby theta nested gamma phase amplitude coupling, by targeting currents active in fast spiking interneurons in the subthreshold regime to tip the balance toward outward currents.
The balance of inward and outward currents in the voltage range traversed between action potentials during the interspike interval (ISI) determines whether a neuron exhibits type 1 or type 2 excitability, and is easily modifiable by altering this balance [55,56]. For type 2, outward currents predominate at steady state in this region of subthreshold membrane potential, but are activated more slowly than the inward currents. Therefore, the dwell time in this region of membrane must be brief so that the outward current does not equilibrate and stop the depolarization caused by the inward currents between spikes. There is a maximum ISI for which repetitive spiking can be supported, hence a minimum frequency below which the neuron cannot sustain repetitive firing.
In contrast, inward currents predominate at steady state in the subthreshold range on membrane potentials spanned by the interspike interval, therefore arbitrarily low firing rates can be sustained.
Decreasing inward or increasing outward currents that have a steady state component during the ISI favors type 2 excitability, and the opposite manipulations favor type 1. These manipulations, depending upon whether inhibition is hyperpolarizing or shunting, could increase theta/gamma PAC. Increased theta/gamma phase amplitude coupling may in turn improve some aspects of cognition.

Network
For all simulations we used 300 neurons of the same excitability type, connected by bi-exponential inhibitory synapses. In the network, the input current for each neuron is given by the following equations:  were uniformly randomly distributed between 0.7 and 3.5 ms.

Numerical simulations and bifurcation analysis
The bifurcation analysis was performed in XPPAUTO [57].
The network models were implemented as Python 2.7 script for the simulation package NEURON [58]. The integration time step was constant at 0.01 ms. Synaptic activation initially was set to zero for all simulations. The membrane potential of each neuron was initialized randomly from a normally distribution with mean -50mV and standard deviation of 20 mV and the slow n variable was initialized at the steady-state value for the membrane potential. The data presented on vector strength, participation, CV of participation and total suppression for networks biased in the oscillatory regime were averaged over ten trials at each parameter setting for runs of duration 2.5 s with the initial 500 ms ignored to minimize the effects of transients. Each trial had its own random connectivity pattern, random initialization of the state variables, random distribution of bias currents, random delay distribution and random noise sources. For the sinusoidal drive simulations the PAC was averaged over 20 periods of sinusoidal drive, again averaged over 10 trials as described above, but no transients were deleted. Simulations were run on Louisiana Optical Network Infrastructure, QB2-cluster. The source code is available at https://senselab.med.yale.edu/Mod-elDB/showmodel.cshtml?model=259366#tabs-1 (will be publicly accessible upon publication; code for reviewers ZhD7rUh64oD1roKUHc7D).
The phase of each spike within a cycle used was calculated using the length of that particular cycle. Cycle lengths are variable and were computed using the peaks in the population rate as described in Tikidji-Hamburyan et al. [17]. The phase was used to construct the vectors for the vector strength calculation.
In the supplementary figures, we applied a multiplicative scale factor ( i F ) to the rates of change of both variables. Since the Gaussian noise term simulates Brownian motion in the membrane potential, in which distance is proportional to the square root of time, the noise term as then divided by the square root of the scale factor. Taken together, these manipulations simply scale the intrinsic frequency since all intrinsic (but not synaptic) processes are sped up or slowed down equally.

Measure of phase amplitude coupling
Since we apply the theta drive, there is no uncertainty with respect to the phase of the theta oscillation, in contrast to the uncertainty in experimental data such as local field potentials or the EEG. This simplified our analysis. Moreover, since we apply the exact same amplitude of theta modulation to different networks, we were interested not only in the tightness of the coupling of the theta drive and the evoked nested gamma oscillation, but also in the magnitude of the gamma oscillations. Therefore in order to quantify the coupling between theta phase and gamma amplitude, we choose the mean vector length (MVL) [59,60], but without normalization of the amplitude.    Table 1, attributes in Table 2) were adjusted to make the comparison as fair as possible. A. Frequency/ current (f/I) curves for type 1 (black dots) and type 2 (red crosses) overlap for the range of bias currents (horizontal bar) used in our heterogeneous networks. The vertical bar shows the range of frequencies exhibited at the heterogeneous bias current levels. The f/I curves were measured using current steps of sufficient duration to allow any transients to die out and establish a steady frequency. Arbitrarily slow frequencies can be obtained for type 1 (not shown) near the bifurcation, but type 2 has a minimum cutoff frequency below which it cannot fire. Between current steps, the membrane potential was returned to its resting value. The bistable range for type 2 is evident from the current values at which a zero frequency quiescent solution coexists with a repetitively firing solution. This region was determined using XPPAUTO. B. The spike shapes and the interspike interval are very similar for type 1 (black curve) and type 2 (red dashed curve).  B1. Minimal inhibition. In contrast to the absence of a fixed point in A1, there is a stable fixed point (technically fixed points, vertical arrow) at the intersection of the average V-nullcline and the n-nullcine. This fixed point is surrounded by an unstable limit cycle (red curve) that forms the boundary between spiking and quiescent trajectories. B2. Inhibition midpoint. As in A2 inhibition accumulates as neurons escape to the right and fire action potentials. The right branch of the Vnullcline is now the firing threshold (vertical arrow). B3. Inhibition maximum. As in A3, two groups are evident corresponding to those that fired and those that skipped. C. Summary. C1. This figure shows that the movement of the fixed point as inhibition waxes and wanes is (leftward) in the fast, horizontal direction for type 1. C2. In contrast, the downward and leftward movement of the fixed point for type 2 has a component in the slow direction, which helps equalize the opportunity of slow neurons to fire as compared to fast neurons.  In contrast to the absence of a fixed point in A1, there is a stable fixed point at the intersection of the average V-nullcline and the n-nullcine. This fixed point is again surrounded by an unstable limit cycle (red curve) that forms the boundary between spiking and quiescent trajectories. B2. Inhibition midpoint. As in A2 inhibition accumulates as neurons escape to the right and fire action potentials. However in this case the ghost of the unstable limit cycle unfurls into a quasithreshold (red curve) and separates spiking and skipping trajectories (vertical arrow labeling the red curve as the firing threshold). B3. Inhibition maximum. As in A3, two groups are evident corresponding to those that spiked and those that skipped (diagonal arrow) . C. Summary. C1. This figure shows that there is less leftward movement of the fixed point as inhibition waxes and wanes in the fast, horizontal direction for type 1 compared to in Fig. 4C1. C2. The fixed point for type 2 tends to trap trajectories inside the unstable limit cycle (red curve, see blowup in inset), which enforces skipping. Top: Representative sparsely firing single neuron type 2 traces with subthreshold oscillations due to network activity. Middle: Raster plots for 300 neurons with faster neurons (based on I app) shown at the top. Bottom: Simulated LFP consisting of summated inhibitory currents throughout the network. C. Heatmaps for 2D parameter space of modulation depth (given in terms of the peak of the sinusoidal conductance waveform) and frequency. Asterisks show parameter values from A and B. C1. Vector strengths for type 1 networks. C2. Vector strengths for type 2 networks. C3. Difference of vector strengths for type 2 and type 1 networks. Hot colors indicate a difference greater than zero. The difference is well above zero in the lower lefthand cornerplot meaning that the vector strength for type 2 networks is higher. Figure 8. Speed of Recruitment of Theta-Nest Gamma in Homogeneous Inhibitory Networks with Hyperpolarizing Versus Shunting Inhibition. A. Homogeneous inhibitory network with hyperpolarizing inhibition showing that type 2 neurons are recruited more quickly into thetanested gamma (leftmost vertical gray bar) and persist longer (rightmost gray bar). A1. Type 1. A2. Type 2. B. Homogeneous inhibitory network with hyperpolarizing inhibition showing that type 1 neurons are recruited more quickly into theta-nested gamma (leftmost vertical gray bar) and persist longer (rightmost gray bar). B1. Type 1. B2. Type 2. Figure 9. Theta Phase Gamma Amplitude Modulation of Heterogeneous Inhibitory Interneuronal Networks with Shunting Inhibition. One example of sinusoidal modulation of these networks, with a 5 Hz sinusoidal modulation in inhibitory conductance shown at the top for the same synaptic conductance and noise standard deviation as in Fig.7. A. Type 1 networks modulation. Top: Representative sparsely firing single neuron type 1 traces with subthreshold oscillations due to network activity. Middle: Raster plots for 300 neurons with faster neurons (based on I app) shown at the top. Bottom: Simulated LFP consisting of summated inhibitory currents throughout the network. B. Type 2 networks. Top: Representative sparsely firing single neuron type 2 traces with subthreshold oscillations due to network activity. Middle:. Raster plots for 300 neurons with faster neurons (based on I app) shown at the top. Bottom: Simulated LFP consisting of summated inhibitory currents throughout the network. C. Heatmaps for 2D parameter space of modulation depth (given in terms of the peak of the sinusoidal conductance waveform)and frequency. Asterisk shows the parameter values from A and B. C1. Vector strengths for type 1 networks. C2. Vector strengths for type 2 networks. C3. Difference of vector strengths for type 2 versus type 1 networks. Cool colors indicate a difference less than zero. The difference is below zero in this plot, meaning that the vector strength for type 2 networks is higher, except at the bottom where the vector strength is relatively small.  Ratio of vector strengths for type 2 versus type 1 networks. Cool colors indicate a ratio less than one. The ratio is always above one in this plot meaning that the vector strength for type 2 networks is higher, as in Figure 7. The additional heterogeneity is the same as in Supplementary  Figure 1.
Supplementary Figure 4. Theta Phase Gamma Amplitude Modulation of Heterogeneous Inhibitory Interneuronal Networks with Hyperpolarizing Inhibition with Additional Heterogeneity in Intrinsic Dynamics. Heatmaps for 2D parameter space of modulation depth and frequency. A. Vector strengths for type 1 networks. B. Vector strengths for type 2 networks. C. Ratio of vector strengths for type 2 versus type 1 networks. Cool colors indicate a ratio less than one. The ratio is always below one in this plot as in Figure 9, meaning that the vector strength for type 2 networks is higher, except at the bottom where the vector strength is so small that there is negligible PAC. The additional heterogeneity is the same as in Supplementary Figure 1.