## Abstract

Fast neuronal oscillations (gamma, 20–80 Hz) have been observed in the neocortex and hippocampus during behavioral arousal. Using computer simulations, we investigated the hypothesis that such rhythmic activity can emerge in a random network of interconnected GABAergic fast-spiking interneurons. Specific conditions for the population synchronization, on properties of single cells and the circuit, were identified. These include the following: (1) that the amplitude of spike afterhyperpolarization be above the GABA_{A} synaptic reversal potential; (2) that the ratio between the synaptic decay time constant and the oscillation period be sufficiently large; (3) that the effects of heterogeneities be modest because of a steep frequency–current relationship of fast-spiking neurons. Furthermore, using a *population coherence* measure, based on coincident firings of neural pairs, it is demonstrated that large-scale network synchronization requires a critical (minimal) average number of synaptic contacts per cell, which is not sensitive to the network size.

By changing the GABA_{A} synaptic maximal conductance, synaptic decay time constant, or the mean external excitatory drive to the network, the neuronal firing frequencies were gradually and monotonically varied. By contrast, the network synchronization was found to be high only within a frequency band coinciding with the gamma (20–80 Hz) range. We conclude that the GABA_{A} synaptic transmission provides a suitable mechanism for synchronized gamma oscillations in a sparsely connected network of fast-spiking interneurons. In turn, the interneuronal network can presumably maintain subthreshold oscillations in principal cell populations and serve to synchronize discharges of spatially distributed neurons.

Although fast gamma cortical oscillation has been the subject of active investigation in recent years (cf. Singer and Gray, 1995), its underlying neuronal mechanisms remain elusive. Two major issues are the cellular origin of rhythmicity (Llinás et al., 1991; McCormick et al., 1993; Wang, 1993) and the mechanism(s) of large-scale population synchronicity (Freeman, 1975; Bush and Douglas 1991; Engel et al., 1991; Hansel and Sompolinsky, 1996). Traditionally, recurrent excitation between principal (pyramidal) neurons is viewed as a major source of rhythmogenesis as well as neuronal synchronization. However, in model studies in which quantitative data about the synaptic time course were incorporated, it was found that glutamatergic synaptic excitation of the AMPA type usually desynchronizes rather than synchronizes repetitive spike firings of mutually coupled neurons (Hansel et al., 1995; van Vreeswijk et al., 1995). Therefore, recurrent connections between pyramidal cells alone do not seem to account for the network coherence during cortical gamma oscillations. It was suggested that pyramidal cell populations may be entrained by synchronous rhythmic inhibition originating from fast-spiking interneurons (Buzsáki et al., 1983; Lytton and Sejnowski, 1991). During field gamma oscillations, intracellular recordings from pyramidal cells revealed both EPSPs and IPSPs phase-locked to the field oscillation frequencies (Jagadeesh et al., 1992; Chen and Fetz, 1993;Soltész and Deschênes, 1993).

In this paper, we address the question whether, in the hippocampus, an interneuronal network can generate a coherent oscillatory output to the pyramidal neurons, thereby providing a substrate for the synaptic organization of coherent gamma population oscillations. In the behaving rat, physiologically identified interneurons were shown to fire spikes in the gamma frequency range and phase-locked to the local field waves (Bragin et al., 1995). Intracellular studies and immunochemical staining demonstrated that these interneurons are interconnected via GABAergic synapses (Lacaille et al., 1987; Sik et al., 1995;Gulyás et al., 1996). Theoretical studies suggest that these GABAergic interconnections may synchronize an interneuronal network when appropriate conditions on the time course of synaptic transmission are satisfied (Wang and Rinzel, 1992, 1993; van Vreeswijk et al., 1995). Moreover, in a recent *in vitro* experiment (Whittington et al., 1995; Traub et al., 1996), the excitatory glutamate AMPA and NMDA synaptic transmissions were blocked in the hippocampal slice. When metabotropic glutamate receptors were activated, transient oscillatory IPSPs in the 40 Hz frequency range were observed in pyramidal cells. These IPSPs were assumed to originate from the firing activities of fast-spiking interneurons synchronized by their interconnections. Computer simulations (Whittington et al., 1995;Traub et al., 1996) lend further support to this hypothesis.

To assess whether an interneuronal network can subserve an adequate basis for the gamma frequency population rhythm in the hippocampus, it is necessary to identify its specific requirements on the cellular properties and network connectivities, as well as to determine whether these conditions are satisfied by particular interneuronal subtypes. The present study is devoted to investigate such requirements using computer simulations. We found that synaptic transmission via GABA_{A} receptors in a sparsely connected network of model interneurons can provide a mechanism for gamma frequency oscillations, and we compared the modeling results with the anatomical and electrophysiological data from hippocampal fast spiking interneurons.

## MATERIALS AND METHODS

*Model neuron.* Each interneuron is described by a single compartment and obeys the current balance equation:
Equation 2where C_{m} = 1 μF/cm^{2} and I_{app} is the injected current (in μA/cm^{2}). The leak current I_{L} = g_{L}(V − E_{L}) has a conductance g_{L} = 0.1 mS/cm^{2}, so that the passive time constant τ_{0} = C_{m}/g_{L} = 10 msec; E_{L} = −65 mV.

The spike-generating Na^{+} and K^{+}voltage-dependent ion currents (I_{Na} and I_{K}) are of the Hodgkin–Huxley type (Hodgkin and Huxley, 1952). The transient sodium current I_{Na} = g_{Na}m^{3}_{∞}h(V − E_{Na}), where the activation variable m is assumed fast and substituted by its steady-state function m_{∞} = α_{m}/(α_{m} + β_{m}); α_{m}(V) = −0.1(V + 35)/(exp(−0.1(V + 35)) − 1), β_{m}(V) = 4exp(−(V + 60)/18). The inactivation variable h obeys a first-order kinetics:
Equation 2where α_{h}(V) = 0.07 exp(−(V + 58)/20) and β_{h}(V) = 1/(exp(−0.1(V + 28)) + 1). g_{Na} = 35 mS/cm^{2}; E_{Na} = 55 mV, φ = 5.

The delayed rectifier I_{K} = g_{K}n^{4} (V − E_{K}), where the activation variable n obeys the following equation:
Equation 2with α_{n}(V) = −0.01(V + 34)/(exp(−0.1(V + 34)) − 1) and β_{n}(V) = 0.125 exp(−(V + 44)/80); g_{K} = 9 mS/cm^{2}, and E_{K} = −90 mV.

These kinetics and maximal conductances are modified from Hodgkin and Huxley (1952), so that our neuron model displays two salient properties of hippocampal and neocortical fast-spiking interneurons. First, the action potential in these cells is followed by a brief afterhyperpolarization (AHP) of about −15 mV measured from the spike threshold of approximately −55 mV (McCormick et al., 1985; Lacaille and Williams, 1990; Morin et al., 1995; Zhang and McBain, 1995). Thus, during the spike repolarization the membrane potential reaches a minimum of about −70 mV, rather than being close to the reversal potential of the K^{+} current, E_{K} = −90 mV. This is accomplished in the model by relatively small maximal conductance g_{K} and fast gating process of I_{K} so that it deactivates quickly during spike repolarization.

Second, these interneurons have the ability to fire repetitive spikes at high frequencies (with the frequency–current slope up to 200–400 Hz/nA) (McCormick et al., 1985; Lacaille and Williams, 1990; Zhang and McBain, 1995). With fast kinetics of the inactivation (h) of I_{Na}, the activation (n) of I_{K}, and the relatively high threshold of I_{K}, the model interneuron displays a large range of repetitive spiking frequencies in response to a constant injected current I_{app} (Fig.1*A*, *left*). It has a small current threshold (the rheobase I_{app} ≃ 0.2 μA/cm^{2}), and the firing rate is as high as 400 Hz for I_{app} ≃ 20 mA/cm^{2}. Similar to cortical interneurons (McCormick et al., 1985; Lacaille and Williams, 1990), the whole frequency–current curve is not linear, and the frequency–current slope is larger at smaller I_{app} values (lower frequencies) (Fig.1*A*, *right*). As a consequence, the neural population is more sensitive to input heterogeneities at smaller I_{app} values. This is demonstrated in Figure1*B*, where a Gaussian distribution of I_{app} is applied to a population of uncoupled neurons (N = 100), with a mean I_{μ} and standard deviation I_{ς}. Given a fixed and small I_{ς} = 0.03, the mean drive I_{μ} is varied systematically, and the resulting dispersion in the neuronal firing frequencies, f_{ς}/f_{μ} (standard deviation of the firing rate/mean firing rate) is shown as function of I_{μ} (Fig. 1*B*,*top*). When plotted versus f_{μ}, it is evident that with the same amount of dispersion in applied current (I_{ς}) the dispersion in firing rates f_{ς}/f_{μ} is dramatically increased for f_{μ} < 20 Hz (Fig.1*B*, *bottom*). This feature has important implications for the frequency-dependent network behaviors (see Results).

*Model synapse.* The synaptic current I_{syn} = g_{syn}s(V − E_{syn}), where g_{syn} is the maximal synaptic conductance and E_{syn} is the reversal potential. Typically, we set g_{syn} = 0.1 mS/cm^{2} and E_{syn} = −75 mV (Buhl et al., 1995). The gating variable s represents the fraction of open synaptic ion channels. We assume that s obeys a first-order kinetics (Perkel et al., 1981; Wang and Rinzel 1993):
Equation 2where the normalized concentration of the postsynaptic transmitter-receptor complex, F(V_{pre}), is assumed to be an instantaneous and sigmoid function of the presynaptic membrane potential, F(V_{pre}) = 1/(1 + exp(−(V_{pre} − θ_{syn})/2)), where θ_{syn} (set to 0 mV) is high enough so that the transmitter release occurs only when the presynaptic cell emits a spike. The channel opening rate α = 12 msec^{−1} assures a fast rise of the I_{syn}, and the channel closing rate β is the inverse of the decay time constant of the I_{syn}; typically, we set β = 0.1 msec^{−1} (τ_{syn} = 10 msec). An example of I_{syn} and IPSP elicited by a single presynaptic spike is illustrated in Figure 1*C*.

*Random network connectivity.* The network model consists of N cells. The coupling between neurons is randomly assigned, with a fixed average number of synaptic inputs per neuron, M_{syn}. The probability that a pair of neurons are connected in either direction is p = M_{syn}/N. For comparison, we also used fully coupled (all-to-all) connectivity (M_{syn} = N). In the model, the maximal synaptic conductance g_{syn} is divided by M_{syn}, so that when the number of synapses M_{syn} is varied, the total synaptic drive per cell in average remains the same.

M_{syn} is the convergence/divergence factor of the neural coupling in the network. Experimentally, an estimate of this important quantity has been obtained for an interneuronal network of the CA1 hippocampus (Sik et al., 1995). A parvalbumin-positive (PV^{+}) basket interneuron was stained intracellularly by biocytin *in vivo*. Its axonal arborization was largely confined in the striatum pyramidale (Fig.2*A*). Other PV^{+}interneurons were stained immunochemically, and the contacts made by the biocytin-filled cell with other PV^{+} cells were counted (Sik et al., 1995). It was concluded that a single PV^{+}basket cell makes synaptic contacts with at least 60 other PV^{+} cells (a majority of which are basket cells) within a spatial region of the volume up to 0.1–0.2 mm^{3} (Sik et al., 1995) (Fig. 2*B*). This tissue volume contains as many as 500–1000 other PV^{+} cells, because the PV^{+} neurons in the pyramidal layer have a cell density of 5.4 × 10^{3} cells/mm^{3} (Aika et al., 1994). Hence, for the CA1 network of basket cells, the experimentally estimated M_{syn} is ∼60. The probability of postsynaptic contacts, however, decreases with the distance between the cell pair (Fig. 2*C*).

*Heterogeneous input.* In the model, single neurons are not identical. Each receives a depolarizing current I_{app} of different intensity and, hence, has a different firing rate. The bias current I_{app} has a Gaussian distribution with a mean I_{μ} and a standard deviation I_{ς}. The parameter I_{μ} determines the mean excitation by the external drive; I_{ς} measures the degree of the heterogeneity in the cell population.

*A measure of network coherence.* To quantify the synchronization of neuronal firings in the network, we introduce a*coherence measure* based on the normalized cross-correlations of neuronal pairs in the network (Gerstein and Kiang, 1960; Welsh et al., 1995). The coherence between two neurons i and j is measured by their cross-correlation of spike trains at zero time lag within a time bin of Δt = τ. More specifically, suppose that a long time interval T is divided into small bins of τ and that two spike trains are given by X(l) = 0 or 1, Y(l) = 0 or 1, l = 1, 2, … , K (T/K = τ). Thus, we define a coherence measure for the pair as:
Equation 2We have also used a slightly different definition where the mean firing rates are substracted from X(l) and Y(l); the substraction did not significantly change our results reported below.

The *population coherence measure* κ(τ) is defined by the average of κ_{ij}(τ) over many pairs of neurons in the network. This coherence measure presents a number of useful properties. First, it is naturally based on cross-correlations and, although we use it here to describe synchronization of oscillations, it can be applied to quantify the synchrony of nonoscillatory neuronal firings. It is calculated from spike trains, requires relatively small sample sizes, and can be used for data analysis of experimental extracellular multiunit recordings. Second, κ(τ) is between 0 and 1 for all τ. For very small τ, κ(τ) is close to 1 (resp. 0) in the case of maximal synchrony (resp. asynchrony). Third, the degree of synchrony of the network dynamics can be quantified by how κ(τ) behaves *as function of* τ. In the case of full synchrony, κ(τ) is 1 for all nonzero τ values; whereas in the case of total asynchrony, κ(τ) is a linearly increasing function of τ (see below). Finally, the *distribution of*κ_{ij}(τ) over neural pairs provides detailed information about the interneuronal interactions and synchronization.

In Results, the population coherence measure κ(τ) is calculated by averaging over all neural pairs in the network of size N. Typically N = 100.

*Numerical methods.* The network model was integrated using the fourth-order Runge–Kutta method, with a time step of 0.05 msec. For random network connectivity and heterogeneity, each set of simulations was run with three to five random realizations of the network connections and applied current distribution. As initial conditions, the membrane potential is uniformly distributed between −70 and −50 mV and the other channel-gating variables are set at their corresponding steady-state values. Coherence was calculated after 1000 msec transients. Simulations were performed on a SUN Sparc Station or a Y-MP Cray Supercomputer.

## RESULTS

### Spike afterhyperpolarization, inhibition, and synchronization

We start by considering a simple case in which all individual cells are identical (i.e., without heterogeneity) and are coupled in an all-to-all fashion (i.e., without randomness in connectivity). As shown in Figure 3*A* (*left*), cells starting at random and asynchronous initial conditions quickly become synchronized and within 5–6 oscillatory cycles their spiking times are perfectly in-phase. The spike AHP amplitude is 15 mV measured from the spike threshold (−52 mV), so the maximal AHP, V_{AHP} = −67 mV, stays above the reversal potential of the synaptic current (E_{syn} = −75 mV). The inequality V_{AHP} > E_{syn} means that during the time course of an action potential and its repolarization, the synaptic action is always hyperpolarizing. This relationship between intrinsic and synaptic properties was found to be an important condition under which the global network synchronization can be achieved (Fig.3*B,C*). In the example (Fig. 3*B*), the speed of the I_{Na} inactivation and the I_{K} activation is slowed down (φ = 3.33 instead of 5), so that repolarization becomes larger (V_{AHP} ≃ −73 mV, close to E_{syn}). In this case, the network synchronization takes much longer time to realize. With φ = 2 (Fig.3*C*), V_{AHP} is −78 mV, which is below E_{syn}, and global network synchrony is lost. Instead, the network is dynamically broken into two clusters: within each cluster the cells fire spikes simultaneously, and the two clusters alternate in time. Such clustering dynamics is a commonly observed behavior of interneuronal networks (Golomb et al., 1994). Hence, synaptic inhibition of the GABA_{A}-type provides a mechanism by which a macroscopic coherence of the network (global synchrony or clustering) can be realized.

To investigate further the dependence of the network synchronization on the inhibitory nature of synaptic interactions, the intrinsic cell properties were kept unchanged, while the reversal potential E_{syn} was gradually varied (Fig.4). In Figure 4*A* the global coherence index κ (compare Eq. 2.5), plotted versus E_{syn}, remains at 1 (perfect synchrony) for E_{syn} < V_{AHP}. It displays an abrupt drop for E_{syn} > V_{AHP} and is close to 0 for E_{syn} > −60 mV. The oscillation frequency dramatically increases as the synaptic effect becomes depolarizing (not shown). Unlike the two cluster dynamics of Figure 3*C* with the coherence index κ = 0.5, in the regime characterized by κ ∼ 0 the relative timing of neuronal firings is essentially arbitrary. This happens in our network model *when the synaptic interactions are excitatory*. In the example given in Figure 4, *B* and *C*, E_{syn} = 0 mV and τ_{syn} = 2 msec, so that the synaptic current mimic that of the glutamate AMPA type. Because recurrent excitation considerably enhances the neural discharge rates, weaker external drive (I_{μ} = 0.1 instead of 1) was used so as to obtain a similar (∼40 Hz) firing frequency with excitatory rather than inhibitory interactions. In this case, although all neurons have very similar rhythmic frequencies, their relative firing phases are almost uniformly distributed between zero and the common oscillation period T (Fig.4*B*). Hence, the probability of coincident firing within a time bin τ between two cells increases proportionally with τ, and the network coherence function κ(τ) grows linearly from zero at τ = 0 to its maximal value of 1 at τ = T msec (Fig. 4*C*).

### Heterogeneity and asynchrony

It is intuitively expected that network synchrony cannot be globally maintained if individual neurons display very different intrinsic oscillation frequencies. We studied the effects of heterogeneity, using a Gaussian distribution of the applied current intensity I_{app} with standard deviation I_{ς}. As illustrated in Figure5*A*, the network coherence deteriorates quite rapidly with increasing I_{ς}. This sensitivity is related to the large frequency–current slope of fast-spiking interneurons (Fig. 1*A,B*), so that in a population of uncoupled cells a small current dispersion may imply a wide distribution of firing frequencies. When neurons are synaptically coupled, the distribution of firing frequencies is a product of their interactions. As shown in Figure 5*B*, for small dispersion in I_{app} (I_{ς} < 0.02), minor differences in intrinsic firing rates are overcome by the coupling, and the standard deviation of firing rates f_{ς} = 0. As the network coherence erodes with larger I_{ς} values, f_{ς}grows linearly with I_{ς} (Fig.5*B*). This linear frequency–current relationship*in standard deviation* is a network property of coupled cells. By contrast, the mean firing rate f_{μ}changes only slightly, illustrating the relative independence of the neural firing rates and synchronicity. The moderate decrease in f_{μ}, however, *is* related to the decreased degree of network coherence, because asynchronized cell firings result in an averaged tonic hyperpolarization that slows down the firing rate (see below).

Figure 5, *C* and *D*, illustrates a partially coherent state. In the rastergram (Fig. 5*C*), neurons are labeled in the increasing order of their I_{app}values. The cells with the smallest injected currents are not in synchrony with those firing at higher rates. The population coherence measure κ(τ) increases sharply with τ and reaches the value of 1 at τ ≃ 1/f_{μ} (Fig. 5*D*). Considered as a function of τ, κ(τ) may be viewed as a distribution function of the neural pairs whose relative firing phase is τ, between 0 and the mean oscillation period (estimated as T_{μ} = 1/f_{μ}). The *derivative*dκ/dτ, therefore, is the corresponding distribution density. Network synchronization is manifested by a peak at zero phase of dκ/dτ (Fig. 5*E*), reflecting the sharp increase of κ near τ = 0. Note that a second peak is expected near τ = T_{μ} (Fig. 5*E*) because the spiking event is periodic.

The network dynamics with I_{ς} = 0.1 is illustrated in Fig. 5,*F– H*. By contrast to Figure5*C–E*, here the rastergram is quite irregular (Fig.5*F*). The linear form of the function κ(τ) (Fig.5*G*) and its flat derivative (Fig.5*H*) are consistent with a totally desynchronized behavior, where the relative firing phase of neural pairs is*uniformly* distributed between 0 and T_{μ}.

Without the aid of the coherence function κ(τ), it would be difficult to conclude from the rastergram (F) that the network is completely disordered. Indeed, if one looks at the summated synaptic drive, s(t) = (1/N) ∑^{N}_{i=1}s_{i}(t), oscillatory fluctuations are significant in this “population field” (Fig. 6*B*,*top*). This is because, if every cell oscillates, the summation of a relatively small number (e.g., N = 100) of oscillatory signals would still show some oscillations, even if cells are completely asynchronous. To assess the macroscopic coherence of the network, one can compute the temporal variance ς^{2} of s(t) for different network sizes and assess whether the fluctuations in s(t) persist in large networks. The network is asynchronous if ς^{2} of s(t) decreases inversely with the network size, ς^{2}(N) ∼ 1/N (Hansel and Sompolinsky, 1996). This is shown to be the case in Figure 6*A*. As N increases, s(t) becomes flatter (Fig.6*B*) and the peak in the power spectrum gets smaller (Fig. 6*C*). This example shows that global network synchrony cannot be assessed correctly by oscillatory fluctuations in the population field by the presence of a peak in the power spectrum if the network size is not sufficiently large.

### Sparse network and minimal connectedness

Random connections among interneurons can be introduced into the model by assuming that a cell makes synaptic contact to a second cell with a probability p. Then, if N is the total number of cells, there are M_{syn} = pN presynaptic cells that converge to a postsynaptic cell, *on average*. As shown above, the network can be synchronized with all-to-all connectivity (M_{syn} = N, p = 1). Because synchronization is impossible without synaptic connections (M_{syn} = 0), we examined the dependence of the population synchrony on M_{syn}. To evaluate the effect of sparse connectivity separately, M_{syn} was gradually varied for a network of 100 identical neurons (i.e., without heterogeneity). As shown in Figure7*A*, the population coherence (as measured by κ) is essentially zero for M_{syn} below a critical value ≃40; above which it starts to become significant, increases rapidly with M_{syn}, and reaches the maximum of 1 in the all-to-all limit (M_{syn} = N). Thus, the dependence of the network synchrony on M_{syn} is highly nonlinear. There is a minimal value of M_{syn} and neural interconnections have to be sufficiently dense to generate global population synchronization. This critical M_{syn} value is not simply a required minimum for the total synaptic drive per cell, because it does not change noticeably when the maximal synaptic conductance is reduced by a factor of 2 (Fig. 7*A*). Also, in the presence of heterogeneity (e.g., with I_{ς} = 0.03), the critical M_{syn} value remains the same, but the quantitative degree of network coherence for M_{syn} larger than the critical value was reduced (data not shown). On the other hand, this minimal connectedness depends on the probability rules in the network design and on the network dynamical state under consideration. For instance, it is much smaller if the number of synapses per cell is exactly the same number M_{syn}, but the actual set of presynaptic cells is chosen randomly (Fig. 7*A*). Or, when the oscillation frequency is increased from ∼35 Hz to ∼100 Hz with I_{μ} = 3 instead of 1, the network is synchronized only with M_{syn} > 75 (Fig.7*A*).

We next examined whether the required minimal connectedness increases in proportion with the network size. Simulations were performed with N = 100, 200, 500, and 1000, and it was found that the onset of network coherence corresponded to a small value of M_{syn}, close to M_{syn} = 60 for large network sizes (Fig. 7*B*). Therefore, the minimal number of synapses per cell that is required for the network synchronization *is not* a fraction of the total number N of cells. It either remains finite for large N or it could conceivably depend weakly on N.

### Partial synchronization in a sparse and heterogeneous network

As stated above, the minimal connectedness for a network coherence remains the same when a moderate amount of heterogeneity is added. Figure 8 illustrates a partially synchronized network, with both sparse connections (M_{syn} = 60) and a dispersion in the external drive (I_{ς} = 0.03). In that case, a major fraction of neurons (group I) displays similar firing rates (close to 39 Hz), although their intrinsic firing rates are different. The remaining neurons (group II) have lower firing rates (below 34 Hz) that are scattered in the diagram. The membrane potential trace of a representive cell is shown in Figure 8*B*, together with its synaptic drive s_{syn} (sum of the synaptic gating variables to that cell) and the population synaptic field, both displaying fairly regular oscillations at the same frequency. Cross-correlations of membrane potentials between the representative cell (a in Fig. 8*C*) and three other cells (*b, c*, and *d* in Fig.8*C*) show that some pairs (ac and ad) are synchronized with near-zero phase shift, but a significant phase difference may be present for other pairs (e.g.,*ab*; Fig. 8*D*).

The normalized cross-correlation of spike trains at zero phase lag, defined as our coherence index κ_{ij} for the pair (i, j), was calculated for all neural pairs in the network and plotted against the difference in the firing rates of the pair ‖f_{i} − f_{j}‖ (Fig.8*E*). Those pairs with both firing rates above 34 Hz (group I) are plotted in the top panel, and the other pairs are plotted in the bottom panel. Comparison of the two panels reveals that high zero-phase synchronization (large κ_{ij}) occurs only in pairs of group I neurons, and with similar firing frequencies (small ‖f_{i} − f_{j}‖). For those pairs in the bottom panel, κ_{ij} is small even for pairs with almost identical firing frequencies. Therefore, the network is subdivided into two groups of neurons, and only group I neurons are synchronized with near zero-phase shift. It is not immediately clear why all neurons in the asynchronous group have lower, but not higher, firing frequencies than the synchronous group.

On the other hand, the neural *pairs* can be classified into three categories, according to whether they are monosynaptically uncoupled, coupled in one direction, or coupled in both directions. Histograms of κ_{ij}, however, do not show conspicious difference between such categories (Fig.8*F*), indicating that the degree of synchronization between a pair is *not* primarily determined by their direct monosynaptic contacts. We conclude, therefore, that interneuronal coherence emerges as a network phenomenon. The global character of network synchrony is quantified by κ, the average of κ_{ij} over all pairs. The population coherence function κ(τ) increases rapidly with the time bin τ and reaches the value of 1 for τ ≃ T_{μ}, where again T_{μ} is the average oscillation period (Fig.8*G*). Its derivative shows a clear peak at τ = 0 (Fig. 8*H*).

The partially synchronous dynamics cannot be maintained if the interneuronal connections are too sparse (Fig. 7). Indeed, when M_{syn} is decreased from 60 to 30, with all of the other parameters remaining the same, the network becomes asynchronous (Fig. 9). In this case, neurons are not locked to a same firing frequency (Fig. 9*A*), the rastergram looks irregular (Fig. 9*C*), the population synaptic field is almost constant in time (although the synaptic drive to a single cell still shows residual oscillatory fluctuations) (Fig.9*B*), and the cross-correlations between cell pairs are flat (Fig. 9*D*). Corroboratively, the zero-phase synchronization is very low for all neural pairs (Fig.9*E,F*), in contrast with Figure 8, *E* and*F*. The global coherence function κ(τ) is linear in τ (Fig. 9*G*), and its derivative is flat (Fig.9*H*), as expected for a fully asynchronous neural network.

### Dependence on the synaptic time constant

It was shown that the 40 Hz oscillations in hippocampal interneurons are sensitive to the decay time constant τ_{syn} of the GABA_{A} synapse (Whittington et al., 1995; Traub et al., 1996). We have confirmed their result that an increased τ_{syn} induces a decrease in the oscillation frequency (Fig. 10*A*,*left*). This frequency reduction occurs largely because the slowly decaying synaptic inhibition accumulates in time, resulting in a *tonic* level of hyperpolarization that counteracts the external depolarization. This is shown in Figure10*B*, in which the synaptic drive to a representative cell is displayed for three different values of τ_{syn}. The time average is indicated by a horizontal line, which is higher (0.5, 0.64, 0.7) for larger τ_{syn} values (10, 20, and 40 msec, respectively).

We also considered how the network synchronization depends on τ_{syn}. Within a fixed time window, a pair of cells has a higher chance to fire simultaneously if their firing frequencies are higher. Hence, to make a meaningful comparison, the coherence measure κ should be calculated with different time bin τ when τ_{syn} is varied. In Figure 10*A*(*right*), we chose to use τ equal to one-tenth the mean oscillation period. It was found that the coherence index displays a peaked region around τ_{syn} = 7 msec; the network synchronization is lost for both smaller and larger τ_{syn}values.

The decrease in the coherence index κ with larger τ_{syn}is a consequence of network heterogeneities. Indeed, as shown in Figure10*A*, if the dispersion in applied current is absent (I_{ς} = 0), or if the connectivity is not random (M_{syn} = 100), κ is larger but still decreases with τ_{syn}. But if both I_{ς} = 0 and M_{syn} = 100, κ (hence, the network synchrony) is now maximal for all larger τ_{syn} values. Note that the heterogeneity in either the applied current or the connectivity produces stronger effects for larger τ_{syn}values, because in the lower frequency range the network is more sensitive to input heterogeneities (see Fig.1*B*).

On the other hand, the decrease in the network synchrony for small τ_{syn} is presumably attributable to a synaptic (dynamical) effect: inhibition should not be too fast compared to the oscillation period, so as to synchronize the network (Wang and Rinzel 1992, 1993;van Vreeswijk et al., 1995; Ermentrout, 1996). This is illustrated in Figure 10*C* with τ_{syn} = 2 msec, where the rastergram is quite irregular (*top*). This is probably related to the fact that with global connectivity (M_{syn} = 100) and without heterogeneity (I_{ς} = 0), the network displays two-clusters (*bottom*, Fig. 10*D*), similar to Figure 3*C*. In that case, the global network synchrony (one-cluster) was also observed with different initial conditions, but it was very sensitive to the network heterogeneity (data not shown).

It follows from the left and right panels of Figure10*A* that a peaked region for κ versus τ_{syn} implies that the network coherence is high only within a limited range of the mean oscillation frequencies. This interesting observation was confirmed with different levels of the external drive I_{μ} (Fig.11*A*). When I_{μ} is larger, the oscillation frequencies are higher (*left*). Moreover, the maximum of κ is located at a smaller value of τ_{syn}(*right*). The network synchrony does not depend simply on τ_{syn}, but on the ratio between τ_{syn} and the mean oscillation period T_{μ}, which is an increasing function of τ_{syn} (Fig. 11*B*). When the coherence index is plotted versus τ_{syn}/T_{μ}, it is small for small τ_{syn}/T_{μ} ratios and peaks at τ_{syn}/T_{μ} ≃ 0.2 for all three external excitation levels (I_{app} = 1, 2, 3; Fig. 11*C*).

### Optimal synchronization in the gamma frequency range

That network synchronization is highest in a limited frequency range appears to be a robust finding in our simulations. For instance, the phenomenon was also observed when the mean input current (I_{μ}) was varied continuously (Fig.12*A*). With stronger external drive, the mean oscillation frequency increases monotonically (Fig. 12*A*, *top*), but the network coherence displays high values only for an intermediate I_{μ} range (Fig. 12*A*,*bottom*). In an all-to-all network of identical neurons, the network coherence was found to remain maximal (κ = 1) for the entire I_{μ} range (data not shown). Furthermore, at small I_{μ}, heterogeneity in an all-to-all network reduces the synchrony dramatically (with I_{μ} = 0.4, κ = 1, and 0.1 for I_{ς} = 0 and 0.03, respectively). Thus, the decrease of the network synchrony at smaller I_{μ} values is attributable to the high network sensitivity to heterogeneities at lower frequencies (see Fig.1*B*) in addition to a decreased τ_{syn}/T_{μ} ratio (dynamical effect). On the other hand, at larger I_{μ} (higher frequencies), the network coherence requires tighter connectedness and our fixed M_{syn} = 60 may no longer be sufficient (see Fig.7*A*).

We also varied the coupling strength g_{syn}systematically, with three different I_{μ}values. With stronger synaptic inhibition, the oscillation frequency decreases monotonically (Fig. 12*B*, *top*). The network coherence, on the other hand, shows a peaked region at intermediate g_{syn} values (Fig.12*B*, *bottom*). Again, without input heterogeneity and coupling sparseness, the network coherence is maximal (κ = 1) for the entire g_{syn} range (data not shown). Therefore, similar to the case of I_{μ}variation, the decrease of κ is presumably caused by a reduced stability of the network synchrony against input heterogeneity at low frequencies (large g_{syn}) and by a lack of sufficiently tight connectivity (combined with a weakened coupling strength) at high frequencies (small g_{syn}).

When the network coherence index is plotted against the mean frequency, for all four curves in Figure 12, *A* and *B*, it is clear that the network synchronization is realized only in a relatively narrow range of oscillation frequencies (30–80 Hz; Fig.12*C*). As shown by the three curves from Figure12*B*, this frequency range of high network synchrony can be shifted by the level of external drive. With I_{μ} = 2 or 3 (instead of 1), the network is more excited and the coherence peak is located at higher frequencies, as expected, and the peak is somewhat enlarged. However, the amplitude of the peak is dramatically reduced (compare I_{μ} = 1 and 3). This can be explained, again, by the fact that at higher frequencies, the network synchronization requires denser connections (see Fig. 7*A*), whereas here M_{syn} = 60 is kept constant. Figure12*C* demonstrates that, although the optimal frequency range for synchronization is not precisely defined and does depend quantitatively on network parameters and external drive, the high network coherence is robustly limited to a frequency band that coincides roughly with the gamma (20–80 Hz) frequency range.

## DISCUSSION

We examined the emergence of synchronous gamma oscillations in an interneuronal network model. The following conditions were identified for the synchronizing mechanism by GABA_{A} synaptic inhibition. (1) The spike afterhyperpolarization of single neurons should be above the synaptic reversal potential, so that the effect of synaptic inputs is always hyperpolarizing. (2) The synaptic current decay should be relatively slow, such that the ratio between the decay time constant and the oscillation period is not small. (3) Heterogeneities should be sufficiently small. (4) The random network connectedness should be higher than a well defined minimum, which is not sensitive to the network size. When the four conditions are fulfilled, a large-scale network coherence was observed preferentially in the gamma (20–80 Hz) frequency range, although uncoupled neurons are potentially capable of discharging in a wide range of frequencies (0–400 Hz).

### Synchronization by synaptic inhibition

Rhythmogenesis in many biological central pattern generators is based on circuits of coupled inhibitory neurons exhibiting postinhibitory rebound excitation (Selverston and Moulins, 1985). In such systems, however, rhythmic patterns are usually slower than the kinetic time constants of the inhibitory synapses, and reciprocally connected neurons typically fire out-of-phase (Perkel and Mulloney, 1974). More recently, it was recognized that neural oscillations can be synchronized by mutual inhibition at zero phase shift, provided that the synaptic kinetics is sufficiently slow as compared to the oscillation period (Wang and Rinzel, 1992, 1993). Intuitively, slow synaptic decay offers the possibility for neurons to “escape” synchronously as the synaptic inhibition wanes below a certain threshold, thus fire at the same time. A general conclusion from this scenario is that synapses with given temporal characteristics may be suitable to synchronizing large neural populations in a particular oscillatory frequency range; a major determining factor is the ratio between the synaptic time constant and the oscillation period. Other computational works have since shown that the mechanism of synchronization by inhibition may be quite general (van Vreeswijk et al., 1995; Whittington et al., 1995; Ermentrout, 1996; Traub et al., 1996). Synchronization by GABA_{A} synapses is facilitated if the synaptic reversal potential is more negative than the maximum spike afterhyperpolarization. This, however, is necessary only for the perfect global, but not for partial, synchronization. Although fast-spiking interneurons typically display larger AHP amplitudes than pyramidal cells, their measured maximal AHP is usually not below −70 mV (McCormick et al., 1985; Lacaille and Williams, 1990; Morin et al., 1995; Sik et al., 1995), hence, above the reversal potential of −75 mV for GABA_{A} synapses (Buhl et al., 1995). On the other hand, the time constant τ_{syn} of the synaptic current determines the range of the oscillation frequencies where the synchronization can be realized by mutual inhibition. For gamma oscillations (∼40 Hz), the τ_{syn} should be larger than 3 msec, compatible with the estimated τ_{syn} (Otis and Mody, 1992). Therefore, these requirements seem to be fulfilled by cortical fast-spiking interneurons, especially the basket cells, interconnected by GABA_{A} synapses. On the other hand, the synchronization mechanism does not require interneurons to be endowed with a postinhibitory rebound property.

The network synchronization was quantified by a coherence index that was defined in terms of the zero-time cross-correlations of spike trains. This proved to be a useful and reliable measure of population synchrony. By contrast, population-averaged quantities like the “synaptic field” may display significant oscillatory fluctuations even when most neurons are in fact asynchronous, *if* the size of the probed neural population is small. In that case, however, the field fluctuations decrease with the network size and become almost flat for large networks. This observation suggests caution in the interpretation of oscillatory local field potential in experimental measurements and in small-network simulations.

We studied the dependence of the network coherence on the heterogeneity in interneuronal properties. Typically, with moderate heterogeneities, coupled oscillators break down into a synchronous and an asynchronous subpopulations. In our model, asynchronous neurons in such a partially synchronous state all have lower firing frequencies than the synchronous ones. The observed high sensitivity on the degree of heterogeneity can be attributed partly to the large frequency–current slope of these fast-spiking cells, so that a minor dispersion in the external drives may result in a wide distribution of single cells’ oscillation frequencies. It would be of interest to investigate whether the network synchronization may become more robust in the presence of heterogeneities and noise, if the neural connections are structured in space (Somers and Kopell, 1995), or when an excitatory pyramidal population is included that reciprocally interacts with the interneuronal population (Kopell and LeMasson, 1994; Wang et al., 1995). Moreover, information processing in the cortex may involve a small subset of pyramidal neurons at a time. These pyramidal cells, by activating their common interneuron targets, may exert localized effects on the synchronous oscillations in a selective neural subpopulation of the cortical network.

### Random connection and critical sparseness

One of our main objectives was to assess how dense synaptic interconnections must be for the maintenance of the synchronized network oscillations. The degree of network coherence depended on the average number of synaptic contacts per cell, M_{syn}, in a highly nonlinear fashion. A minimal M_{syn} value was identified, below which the network becomes totally asynchronous. For M_{syn}above its critical value, the degree of network coherence becomes nonzero and increases with M_{syn}. We demonstrated that this critical connectedness does not increase in proportion with the number of neurons, hence, is not sensitive to the network size. However, further analysis is needed to determine whether there is a very weak (e.g., logarithmic) dependence. Such a dependence is expected, for instance, if the required connectedness is related to the minimal number of links for a random network to be topologically connected at large scales (Erdös and Rényi, 1960). The existence of a small critical connectedness has also been found in other neural network models (Barkai et al., 1990; Wang et al., 1995), suggesting that this may be a general feature of sparsely connected random neural networks.

Note that the actual value of this minimal connectedness depends on the details of single-cell and network properties, as well as on the cooperative dynamics under consideration. Indeed, we found that synchronization of oscillations at higher frequencies requires tighter interneuronal connections (Fig. 7*A*). Further analyses are needed to provide a theoretical understanding of this simulation result. For rhythmicities in the gamma frequency range, our model network require 60 or more synaptic contacts from an interneuron to other interneurons. This number is comparable with the estimated divergence/convergence factor of CA1 basket cells in the hippocampus (Sik et al., 1995).

We also demonstrated that synchronization of oscillatory neurons usually cannot be attributed simply to the presence of monosynaptic couplings between the respective cells. The macroscopic network synchronization is thus a truly emergent phenomenon of large neuronal aggregates.

### Frequency selection for synchronization

An important finding of the present study is that interneuronal networks can be synchronized by GABA_{A} synapses preferentially within the gamma frequency range. This happens despite the wide range (0–400 Hz) of possible single-neuron firing rates. For instance, when the synaptic time constant (τ_{syn}), external drive (I_{μ}), or coupling strength (g_{syn}) is varied gradually, the neuronal firing frequencies change monotonically and cover a wide frequency range. In all three cases, however, the degree of network synchronization shows a relatively narrow peak within the gamma frequency range (20–80 Hz).

This phenomenon of frequency selection may be qualitatively understood in terms of the following three neural and network properties: (1) the high network sensitivity to heterogeneities at low frequencies, attributable to the steep frequency/current curve of single neurons; (2) the minimal connectedness for the synchrony, which is larger at higher oscillation frequencies; (3) the “dynamical effect,” namely, network synchronization is impossible or fragile against heterogeneities if the ratio between τ_{syn} and the oscillation period T (τ_{syn}/T) becomes too small. Hence, given τ_{syn} ≃ 10 msec, if the oscillation frequency is lower than 20 Hz (T > 50 msec), the network coherence may be abolished by a high sensitivity to input heterogeneity (see Fig. 1*B*) and by a reduced τ_{syn}/T ratio. On the other hand, if the frequency is higher than 80 Hz, the synchrony may no longer be possible because the interneuronal connectivity is not sufficiently tight (see Fig. 7*A*). Consequently, the network can be highly synchronized only at 20–80 Hz.

We would like to emphasize that this frequency band for coherent oscillations cannot be determined in an absolute precision, and its quantitative limits do depend on the details of the model. However, our results robustly demonstrated that the synaptic time constant delimits a frequency band of coherent network oscillations, and the GABA_{A}-type synapse (with τ_{syn} ≃ 10 msec) seems especially suitable for the gamma rhythmicity. The limited frequency range of population gamma oscillations has been observed in the behaving rat (Bragin et al., 1995) and in hippocampal slices (Whittington et al., 1995; Traub et al., 1996). Our findings suggest that the synchronization mechanism by an interneuronal network would be especially effective if the fast oscillations are generated in the gamma frequency range by pacemaker neurons (Llinás et al., 1991;McCormick et al., 1993). On the other hand, even in the absence of pacemaker neurons, coherent field oscillations could still be observed in the gamma frequency range, not because possible firing rates of inhibitory interneurons are restricted to a narrow frequency band, but because outside this frequency range the large-scale coherence is not possible (by this particular synaptic mechanism).

### Physiological implications

Recurrent excitatory connections have long been regarded as a possible synaptic substrate underlying correlated neural firings in general and massively synchronous brain rhythms in particular. Although coherent slow (epileptic) rhythmic bursting may indeed emerge in a disinhibited pyramidal cell network (Chagnac-Amitai and Connors, 1989;Traub et al., 1993), modeling studies suggest that mutual excitation via the AMPA-type synapse often cannot synchronize neural oscillations in the gamma frequency range, at least for simple repetitive spiking neurons (Hansel et al., 1995; van Vreeswijk, et al., 1995) (present work). An alternative mechanism for the fast entrainment of principal cells has emerged recently. Previous work (Bragin et al., 1995;Whittington et al., 1995; Traub et al., 1996) and the present model suggest that networks of interneurons are critically involved in the generation of coherent gamma oscillations. The advantage of such “synchronizer” function of interneuronal networks is the maintenance of subthreshold and coherent modulation of the large, sparsely connected pyramidal cell populations and the resulting precise timing of their action potentials (Buzsáki and Chrobak, 1995;Hopfield, 1995).

## Footnotes

This work was supported by the National Institute of Mental Health (MH53717-01), Office of Naval Research (N00014-95-1-0319), and the Sloan Foundation to X.J.W.; and HFSP and the National Institute of Neurological Disease and Stroke (NS34994) to G.B. and X.J.W. Simulations were partly performed at the Pittsburgh Supercomputing Center. We thank D. Golomb, D. Hansel, J.-C. Lacaille, and C. McBain for discussions, A. Sik for preparing Figure 2, and L. Abbott, J. Lisman, and R. Traub for carefully reading this manuscript.

Correspondence should be addressed to Xiao-Jing Wang, Center for Complex Systems, Brandeis University, Waltham, MA 02254.