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 GABAA 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 GABAA 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 GABAA 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 GABAA 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 Cm = 1 μF/cm2 and Iapp is the injected current (in μA/cm2). The leak current IL = gL(V − EL) has a conductance gL = 0.1 mS/cm2, so that the passive time constant τ0 = Cm/gL = 10 msec; EL = −65 mV.
The spike-generating Na+ and K+voltage-dependent ion currents (INa and IK) are of the Hodgkin–Huxley type (Hodgkin and Huxley, 1952). The transient sodium current INa = gNam3 ∞h(V − ENa), 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). gNa = 35 mS/cm2; ENa = 55 mV, φ = 5.
The delayed rectifier IK = gKn4 (V − EK), 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); gK = 9 mS/cm2, and EK = −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, EK = −90 mV. This is accomplished in the model by relatively small maximal conductance gK and fast gating process of IK 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 INa, the activation (n) of IK, and the relatively high threshold of IK, the model interneuron displays a large range of repetitive spiking frequencies in response to a constant injected current Iapp (Fig.1 A, left). It has a small current threshold (the rheobase Iapp ≃ 0.2 μA/cm2), and the firing rate is as high as 400 Hz for Iapp ≃ 20 mA/cm2. 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 Iapp values (lower frequencies) (Fig.1 A, right). As a consequence, the neural population is more sensitive to input heterogeneities at smaller Iapp values. This is demonstrated in Figure1 B, where a Gaussian distribution of Iapp 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 Isyn = gsyns(V − Esyn), where gsyn is the maximal synaptic conductance and Esyn is the reversal potential. Typically, we set gsyn = 0.1 mS/cm2 and Esyn = −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(Vpre), is assumed to be an instantaneous and sigmoid function of the presynaptic membrane potential, F(Vpre) = 1/(1 + exp(−(Vpre − θ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 Isyn, and the channel closing rate β is the inverse of the decay time constant of the Isyn; typically, we set β = 0.1 msec−1 (τsyn = 10 msec). An example of Isyn 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, Msyn. The probability that a pair of neurons are connected in either direction is p = Msyn/N. For comparison, we also used fully coupled (all-to-all) connectivity (Msyn = N). In the model, the maximal synaptic conductance gsyn is divided by Msyn, so that when the number of synapses Msyn is varied, the total synaptic drive per cell in average remains the same.
Msyn 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 mm3 (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 × 103 cells/mm3 (Aika et al., 1994). Hence, for the CA1 network of basket cells, the experimentally estimated Msyn 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 Iapp of different intensity and, hence, has a different firing rate. The bias current Iapp 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 acoherence 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.
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, VAHP = −67 mV, stays above the reversal potential of the synaptic current (Esyn = −75 mV). The inequality VAHP > Esyn 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 INa inactivation and the IK activation is slowed down (φ = 3.33 instead of 5), so that repolarization becomes larger (VAHP ≃ −73 mV, close to Esyn). In this case, the network synchronization takes much longer time to realize. With φ = 2 (Fig.3 C), VAHP is −78 mV, which is below Esyn, 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 GABAA-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 Esyn was gradually varied (Fig.4). In Figure 4 A the global coherence index κ (compare Eq. 2.5), plotted versus Esyn, remains at 1 (perfect synchrony) for Esyn < VAHP. It displays an abrupt drop for Esyn > VAHP and is close to 0 for Esyn > −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, Esyn = 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 Iapp 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 Iapp (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 relationshipin 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 Iappvalues. 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 derivativedκ/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 isuniformly 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=1si(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 Msyn = pN presynaptic cells that converge to a postsynaptic cell, on average. As shown above, the network can be synchronized with all-to-all connectivity (Msyn = N, p = 1). Because synchronization is impossible without synaptic connections (Msyn = 0), we examined the dependence of the population synchrony on Msyn. To evaluate the effect of sparse connectivity separately, Msyn 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 Msyn below a critical value ≃40; above which it starts to become significant, increases rapidly with Msyn, and reaches the maximum of 1 in the all-to-all limit (Msyn = N). Thus, the dependence of the network synchrony on Msyn is highly nonlinear. There is a minimal value of Msyn and neural interconnections have to be sufficiently dense to generate global population synchronization. This critical Msyn 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 Msyn value remains the same, but the quantitative degree of network coherence for Msyn 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 Msyn, 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 Msyn > 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 Msyn, close to Msyn = 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 (Msyn = 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 ssyn (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 ‖fi − fj‖ (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 ‖fi − fj‖). 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 Msyn 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 andF. 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 GABAA 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 τsynvalues.
The decrease in the coherence index κ with larger τsynis 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 (Msyn = 100), κ is larger but still decreases with τsyn. But if both Iς = 0 and Msyn = 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 τsynvalues, 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 (Msyn = 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 (Iapp = 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 Msyn = 60 may no longer be sufficient (see Fig.7 A).
We also varied the coupling strength gsynsystematically, 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 gsyn values (Fig.12 B, bottom). Again, without input heterogeneity and coupling sparseness, the network coherence is maximal (κ = 1) for the entire gsyn 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 gsyn) and by a lack of sufficiently tight connectivity (combined with a weakened coupling strength) at high frequencies (small gsyn).
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 Msyn = 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.
We examined the emergence of synchronous gamma oscillations in an interneuronal network model. The following conditions were identified for the synchronizing mechanism by GABAA 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 GABAA 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 GABAA 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 GABAA 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, Msyn, in a highly nonlinear fashion. A minimal Msyn value was identified, below which the network becomes totally asynchronous. For Msynabove its critical value, the degree of network coherence becomes nonzero and increases with Msyn. 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 GABAA 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 (gsyn) 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 GABAA-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).
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).
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.