Brain rhythms arise through the synchronization of neurons and their entrainment in a regular firing pattern. In this process, networks of reciprocally connected inhibitory neurons are often involved, but what mechanism determines the oscillation frequency is not completely understood. Analytical studies predict that the emerging frequency band is primarily constrained by the decay rate of the unitary IPSC. We observed a new phenomenon of resonant synchronization in computer-simulated networks of inhibitory neurons in which the synaptic current has a delayed onset, reflecting finite spike propagation and synaptic transmission times. At the resonant level of network excitation, all neurons fire synchronously and rhythmically with a period approximately four times the mean delay of the onset of the inhibitory synaptic current. The amplitude and decay time constant of the synaptic current have relatively minor effects on the emerging frequency band. By varying the axonal delay of the inhibitory connections, networks with a realistic synaptic kinetics can be tuned to frequencies from 40 to >200 Hz. This resonance phenomenon arises in heterogeneous networks with, on average, as few as five connections per neuron. We conclude that the delay of the synaptic current is the primary parameter controlling the oscillation frequency of inhibitory networks and propose that delay-induced synchronization is a mechanism for fast brain rhythms that depend on intact inhibitory synaptic transmission.
Networks of inhibitory neurons with reciprocal synaptic connections are found from invertebrate retina (Hartline and Ratliff, 1972) up to mammalian hippocampus (Cobb et al., 1997), striatum, thalamus, cerebellum, and neocortex (Kisvarday et al., 1993; Tamás et al., 1998; Galarreta and Hestrin, 1999; Gibson et al., 1999; McBain and Fisahn, 2001). In many of these brain regions, fast-oscillating local field potentials (from 40 to >200 Hz) were recorded with subclasses of inhibitory neurons firing phase-locked to, and often thought to generate, the rhythm (Bragin et al., 1995; Whittington et al., 1995; Kandel and Buzsáki, 1997; Csicsvari et al., 1999, 2003; Grenier et al., 2001). The rhythms are typically confined within a frequency band characteristic of the brain area, the experimental procedure, or the associated behavior (for review, see Gray, 1994; Buzsáki and Chrobak, 1995; Baçar, 1998; Farmer, 1998; Traub et al., 1999).
If networks of inhibitory neurons are responsible for maintaining the rhythms, then the constancy in frequency in each particular case requires either that the individual neurons have an intrinsic resonance around the oscillation frequency or that populations of inhibitory neurons preferentially synchronize within the recorded frequency band.
Theoretical studies demonstrated that a homogeneous pair of reciprocally connected inhibitory neurons lacking intrinsic resonance can synchronize their spikes with zero phase lag, provided a finite delay separates the generation of an action potential by one of the neurons and the peak of the synaptic response, or IPSC, in the paired neuron (van Vreeswijk et al., 1994; Ernst et al., 1995). Most studies of inhibitory networks, analyzing the parameter dependency of synchrony, implemented the required delay as a slow rise time of the IPSC, thereby avoiding the analytical difficulties imposed by the discrete delays of axonal conduction and synaptic transmission. The frequency selectivity of the rhythms, however, is not completely understood and was attributed to the decay time constant of the synaptic response (Wang and Rinzel, 1993; Traub et al., 1996; Wang and Buzsáki, 1996).
Here we demonstrate that the discrete delay in the onset of the IPSC, which is attributable to the finite duration of axonal conduction and synaptic transmission, induces a resonance phenomenon in inhibitory networks, and the oscillation period at resonance is close to four times the combined axonal and synaptic delay. A delayed IPSC onset also enables synchrony to develop in networks with a fast-rising IPSC waveform, in accordance with recent experimental findings (Bartos et al., 2001, 2002; Carter and Regehr, 2002).
The present study elaborates on an inhibitory network model developed for reproducing 160 Hz oscillations recorded in the cerebellar cortex of transgenic mice deficient for calretinin and calbindin (Cheron et al., 2001; Maex et al., 2002). The observed frequency constancy made us search for a mechanism capable of inducing resonance in spatially organized, heterogeneous inhibitory networks.
Materials and Methods
Model neurons. Fast-spiking model neurons had a spherical soma and an unbranched dendrite of five cylindrical compartments. All compartments had the same surface area and a passive membrane time constant of 30 msec. The spike-generating inactivating Na+ channel and the delayed rectifier were restricted to the soma. All compartments had a high-voltage-activated Ca2+ channel, A-type and combined voltage- and Ca2+-activated K+ channels, and a weak anomalous inward rectifier. Rate constants were converted to 37°C, using Q10 = 3 for voltage-gated and Q10 = 2 for ligand-gated channels (Maex and De Schutter, 1998a). The reversal potential of the leak current was drawn from a uniform distribution between -70 and -60 mV. The neurons did not exhibit subthreshold oscillations or intrinsic membrane resonance, nor did they fire postinhibitory rebound spikes. Suprathreshold current injection evoked narrow spikes (0.3 msec width at -40 mV) without firing rate adaptation (Maex and De Schutter, 1998a, their Fig. 2A). We selected this neuron model, which combines the active properties of a model cerebellar granule cell and the passive properties of a Golgi cell, because its firing rate increased linearly with the intensity of applied current over almost the entire dynamic range (from 15 Hz at the brisk 65 pA threshold to 1000 Hz, the maximum firing rate imposed by a 1 msec absolute refractory period; mean slope, 0.65 Hz/pA for an input resistance of 162 MΩ).
Model GABAA receptor synapses. A GABAA receptor channel with a reversal potential of -70 mV was inserted on the somatic compartment. Each afferent action potential triggered a unitary conductance increase with a fixed, dual-exponential time course: 1 The parameters w, τ, and d are the relative peak conductance, decay time constant and delay of the synaptic response, respectively (see Fig. 2A), and t denotes time elapsed after the rising phase of the afferent action potential crossed the -20 mV level. The parameter d encompasses both axonal conduction and synaptic transmission delays. The normalization constants A and n (the number of afferent synapses) ensured that the summed unitary conductances had a peak amplitude equal to wG when τ or n was varied. The dimensionless parameter w is used only for reference. For w = 1, the summed peak conductance wG was 3 mS/cm2. The rise time constant τr determines the shape of the conductance response and was taken as τ/27.4 to reproduce the IPSC waveform of hippocampal pyramidal neurons (Ropert et al., 1990). For τ = 3 msec, τr = 0.11 msec, which is comparable after temperature correction to the mean rise time constant of 0.16 msec, recorded in basket cells of dentate gyrus at 34°C (Bartos et al., 2002). In selected simulations (see Fig. 7B), neurons reciprocally connected through GABAA receptor synapses also made electrical synapses on the same compartments (Tamás et al., 2000). These gap junctions were pure resistors with a conductance expressed as a fraction of the peak conductance of the associated chemical synapse.
Conversion to physiological data. Because conductances are expressed as densities over the compartments to which the channels are attached, the output of this neuron model is primarily invariant over spatial scaling. For a neuron with total membrane area S, the unitary peak conductance gunitary can be derived as: 2 The factor 6 in the denominator converts density over the somatic compartment to density over the entire model neuron. Hence, taking wG = 3 mS/cm2, the afferent synapses of a neuron with total area S = 12,000 μm2 (Bartos et al., 2001) and n = 60 (Wang and Buzsáki, 1996) have a unitary peak conductance of 1 nS.
Afferent fibers. The dendritic compartments received AMPA receptor synapses (rise and decay time constants, 0.03 and 0.5 msec, respectively; summed peak conductance, G = 16.1 mS/cm2; Maex and De Schutter, 1998a) from a population of afferent fibers. These fibers were Poisson processes, all firing at the same constant rate throughout a simulation. The time histogram of spike counts of the fiber population had a flat power spectrum. Because synapses and fibers are computationally expensive, a diluted population of fibers was simulated. Each fiber represented the compound spike trains from several tens of excitatory fibers and fired at a rate varied between simulations from 50 to 12,800 spikes/sec.
Model network configurations. We first examined how the parameters of the inhibitory synapse affect the oscillation frequency and power in a network characterized by single values for d, τ, and w. We then proceeded to validate the derived relationships on networks with distributed delays and weights and with stochastic connections.
In the one-dimensional network with only nearest neighbor coupling, 100 model neurons were positioned on a one-dimensional array. The neurons made autapses, with zero delay, and synapses on the closest neighbor on each side with delay d. This overtly simplified network has the advantage of being completely characterized by the single delay value d, the single decay time constant τ, and the single synaptic weight value w. We varied systematically the values of d, τ, and w (see Figs. 1, 2, 3A, 7).
In the one-dimensional network with distributed delays, each neuron was connected to all neighbors within an axon or connection radius r (measured in units of interneuron distance) with the delays d set proportional to the interneuron distances. Assuming a spacing of 300 μm between the neurons and an axonal conduction speed of 0.3 m/sec, the delay to the nearest neighbor measured 1 msec, and delays to more distant neurons were an integer multiple of this. The mean delay of the network connections measured (r + 1)/2. We varied systematically the values of r and τ (Figs. 3B, 4).
The two-dimensional network, representing circuits with a predominantly radial connection pattern, was an array of 30 × 30 neurons. Their somata were positioned on a triagonal grid so that each neuron had six equidistant neighbors. Connecting a neuron to all neighbors alike within a radius r makes the delay distribution very skewed and biased toward the largest delay values. Either the connection probability or the synaptic weight, therefore, needed to fall off with distance from the source neuron. In the network documented in Figure 5, synaptic weight decreased exponentially with a space constant twice the interneuron distance. Using r = 8 and a constant connection probability of 0.6, each neuron connected on average to 115 neighbors. The mean delay of the network was the mean d over all connections, each d weighted by the corresponding synaptic weight w and with the exclusion of autapses, which did not contribute to synchrony. We used τ = 3 msec.
For the effects of network size and connectivity (Fig. 6), nontopographic networks were simulated in which all neurons were stochastically connected with a variable probability. All synaptic connections had the same fixed delay, d = 1 msec and τ = 3 msec.
A population of afferent fibers was evenly distributed over the entire network space. The fibers had a conduction speed of 0.3 m/sec and radiated 2.5 mm in a plane perpendicular to the dendritic shafts, making en passant excitatory synapses on the dendritic compartments with a probability of 0.5. There were 8100 fibers in the two-dimensional network and 4000 fibers in the one-dimensional and nontopographic networks (except when stated otherwise), providing each neuron on average with 424 and 320 synapses, respectively.
Boundary conditions. As in our previous studies, normalizing the unitary synaptic conductance over the number of afferent synapses prevented the possibility that the boundary neurons would fire at rates different from the network average. This normalization is justified because it is shown that our main findings are robust to great variability among the neurons in the number of synapses (Fig. 6A). In addition to having fewer synapses, with consequently stronger weights, neurons positioned close to the boundaries had a selective deficit of long connections. The resulting decrease in the mean synaptic delay, which is largest for a neuron positioned at a distance half the connection radius from the boundary, was relatively small. For example, in the one-dimensional network with r = 16, the mean delay of the inhibitory synapses onto neurons 9 and 92 was 7.2 instead of 8.5 msec. Although these inhomogeneities at the boundaries did have a desynchronizing effect on the entire network, the value of the resonance frequency was hardly affected (42.65 Hz in a 100-neuron array, 42.18 Hz for the central 100-neuron segment of a 200-neuron array, and 42.57 Hz in a 100-neuron ring).
Noise and heterogeneity. The networks were noisy and heterogeneous, owing to the randomness of excitation and the randomization of the resting membrane potentials (see above) and synaptic weights (Maex and De Schutter, 1998a). This heterogeneity made the neurons in a disconnected network fire irregularly and with slightly different rates [coefficient of variation (CV) of their firing rate, averaged over all excitation levels, 0.046; CV of their interspike intervals (ISIs), averaged over all neurons, 0.155]. These values typically increased in connected networks because lateral inhibition tended to amplify differences in the firing rate (CV of the firing rate in the one-dimensional network with only nearest neighbor coupling, 0.072; CV of ISIs, 0.210). Finally, a simulation run started from random membrane potentials distributed uniformly between -90 and -20 mV.
Implementation and simulation. Neuronal activity was calculated numerically with a modified version of the Genesis simulator (http://www.genesis-sim.org/GENESIS), using Crank-Nicolson integration in 20 μsec steps. We used a random-number generator with two seed variables to avoid correlations in firing among the afferent fibers.
Analysis of network activity. We monitored network activity as the time series of spike counts in bins of 0.5 msec width (Fig. 1A,D). The power spectrum or periodogram was estimated, using the fast Fourier transform algorithm, as the average of 256-point overlapping Hann (cosine bell) windows (Press et al., 1988). In each window, the mean number of spikes was calculated first and subtracted from the individual counts (or, equivalently, power at f = 0 was set at 0) so that the total power almost equaled the variance of the time signal. Coherent oscillations were quantified by the peak power and the associated center or network frequency. Network frequency was assessed at higher resolution (512 or 1024 point windows) for networks with a low resonance frequency.
To derive the resonance frequency of a network, the network was excited to various levels by incrementing the mean firing rate of the afferent fibers by a factor of the square root of 2 at each new run. The obtained average firing rates of the neurons encompassed their entire dynamic range (10 to >900 spikes/sec). A simulation run produced at least 1000 spikes for each neuron (500 in two-dimensional networks) and had a minimum duration of 3 sec. Over such long runs, power spectra were well reproducible and independent of initial values. From the power spectra obtained at different levels of excitation, a tuning curve was constructed after application of the following normalization procedure. We first noticed that a disconnected network produced a power spectrum with a single peak that was located at the mean neuronal firing rate. This peak had a height almost proportional to the firing rate. In accordance with this, the variance of the time series of spike counts scaled almost linearly with the mean spike count. This relationship between power and firing rate is compatible with a Poisson distribution of spike counts in the bins of the time signal. We therefore divided the peak power obtained at each level of excitation by the associated mean neuronal firing rate, correcting in this way for the power already present in the disconnected network. Finally, peak power scales in a synchronous network quadratically with the number of neurons N and was therefore normalized to a network of N = 100. All frequency-tuning curves plot peak power, divided by the mean neuronal firing rate and by (N/100)2 and multiplied by a scaling factor 1000.
Resonance frequency (fR) is the mean frequency of the tuning curve, which was constructed from at least eight levels of excitation. Taking this weighted average reduces errors attributable to sampling of the tuning curve. If the central peak of the curve was positioned very asymmetrically, the low-power tail was cut to avoid a systematic overestimation and underestimation of fR in networks with low and high fR, respectively. Some tuning curves had a satellite peak at lower frequencies, and fR may have been slightly underestimated in these networks (e.g., see Fig. 6B). Overall, errors on fR are judged to be small because of the systematic attraction of network frequency toward fR (see Fig. 5B).
Peak power is a combined metric of synchrony and rhythmicity and quantifies the network oscillations in a manner comparable with the analysis of experimentally recorded local-field potentials. A metric frequently used in the computational literature is the coherence index (Wang and Buzsáki, 1996; Bartos et al., 2001, 2002) (variations in White et al., 1998; Tiesinga et al., 2000). The coherence index (CI), taking values between 0 and 1, measures exclusively synchrony and assesses, averaged over all neuron pairs of the network, the fraction of coincident spikes. The probability of coincidence increases with the bin width of the discretized spike trains, which by convention is taken in 1/10 of the oscillation period. Measuring only the CI is justified in the above studies because the neurons being excited by constant-current injection were almost perfect oscillators, and the inhibitory synapses constituted the only synchronizing mechanism. In the present study, in contrast, excitation is provided by afferent fibers, which carry random spike trains that can contribute to the synchronization of common efferent neurons (Maex et al., 2000). We mention for each spike raster plot the CI in the figure legend.
We investigated the critical parameters for synchronization and frequency control in a heterogeneous network of model inhibitory neurons, themselves lacking intrinsic resonance. Randomly firing afferent fibers excited the neurons by continually activating AMPA receptor synapses on the dendrites. This noisy but realistic way of driving the neurons accentuated network-induced rhythmicity. Each neuron inhibited its neighbors, within a varying axon radius, through GABAA receptor synapses located on the somatic compartment. The delay d, strength w, and decay time constant τ of the IPSC (see Fig. 2A) were varied, and their effects on network dynamics were compared at various levels of network excitation.
High-frequency oscillations emerge in networks with axonal and synaptic delays
For a better understanding of the resonance mechanism, we first illustrate the results of a one-dimensional array in which a neuron is connected on each side only to the closest neighbor.
Figure 1B demonstrates waxing and waning high-frequency oscillations in a network in which the combined synaptic and conduction delay d was set to 1 msec, which is a typical latency value for inhibitory synaptic currents evoked from nearby neurons in paired recordings (Bartos et al., 2001). As a control, no oscillations were observed in the same network if the onset of the IPSC was instantaneous (delay d = 0; Fig. 1A,C). This latter finding is in accordance with analytical studies proving that synchronous firing is unstable if the rise time of the synaptic response is shorter than the duration of the afferent spike (van Vreeswijk et al., 1994). However, using a synaptic response function with slower rise time compared with this control, e.g., an α-function with instantaneous onset and peak conductance at 1 msec, did not improve synchrony in this sparsely connected network without axonal and synaptic delays.
The oscillation frequency and power were measured at various levels of network excitation spanning the entire dynamic range of the neurons. Figure 1F shows that the tuning curve of the delayed network was centered at a resonance frequency fR of 210 Hz, with power falling off steeply at lower and higher levels of excitation. In the network without delays, no robust oscillations could be evoked at any frequency (Fig. 1E).
At the resonant level of network excitation, the mean firing rate of the neurons was close to the oscillation frequency of the network because each neuron fired a single spike at almost every cycle (Fig. 1D). At nonoptimal levels of excitation, the neurons changed their firing pattern to maintain a network frequency deviating from their firing rate but closer to fR. In particular, at low excitation levels, neurons skipped cycles of the oscillation. As a consequence, adjacent neurons could fire in alternate order (antiphase synchrony) at a rate half the network frequency (Fig. 1F, open data points). At high excitation levels, neurons fired multiple spikes at each single cycle, without affecting much the cycle duration. This doublet or burst firing was most prominent in networks with delays larger than the present 1 msec value (data not shown).
The delay parameter determines the resonance frequency
Several studies of inhibitory networks demonstrated that the oscillation frequency decreases when the strength, decay time constant, or delay of the inhibitory synaptic response are increased (Bush and Sejnowski, 1996; Traub et al., 1996; Wang and Buzsáki, 1996; Pauluis et al., 1999; Bartos et al., 2002; Liley et al., 2002). Figure 2 illustrates the complementary effects of these parameters, i.e., the effects on the frequency-tuning curve of a network.
The frequency tuning was almost invariant over the value w of the synaptic strength (Fig. 2B). In the network with a resonance frequency fR of 210 Hz, illustrated above, a variation in peak conductance by a factor of 8 suppressed the mean firing rate on average by 51%, but fR decreased only 14% (from 225 to 183 Hz). We disregard this weight parameter further and present the tuning curve with largest power for networks simulated at various weight values. Increasing or decreasing the exponential decay time constant τ of the postsynaptic current displaced the tuning curve to lower and higher frequencies, respectively (Fig. 2C). The resonance frequency fR decreased from 338 Hz at τ = 0.75 msec to 149 Hz at τ = 12 msec. This effect was small, however, compared with the inverse relationship obtained between fR and the delay d of the synaptic response (Fig. 2D). Varying d from 0.25 to 4 msec decreased fR from 536 to 65 Hz.
Over a broad range of values of the delay d and the decay time constant τ, fR approximated 1/(4d) (Fig. 3A). At the lowest delay values (0.25 and 0.5 msec), the measured fR was <1/(4d) as the network frequency approached 1000 Hz, the maximal neuronal firing rate. Networks with a (single) large delay value, on the other hand, could have multimodal tuning curves because at high excitation levels, the neurons started firing bursts, in particular when τ was much less than 4d (the oscillation period at resonance). Such unrealistic networks did not always have a clear optimal frequency, and their data were not included in Figure 3A (but see next section).
A different value of τ was optimal for different values of d (Fig. 3A), a nonoptimal τ driving fR off the predicted frequency of 1/(4d) with a considerable loss of power (Fig. 2C). For example, oscillations of ∼250 Hz (obtained with d = 1 msec) were optimal with τ = 0.75-3 msec. Oscillations of ∼62.5 Hz (d = 4 msec) required a τ of 6-12 msec, a value comparable with the recorded 10 msec decay of (probably compound) IPSCs during hippocampal gamma oscillations (Traub et al., 1996). With τ ≤3 msec, networks with only nearest neighbor coupling were unable to generate oscillations <100 Hz.
Distributed delays increase the robustness of frequency control
Increasing the connection radius r of the one-dimensional network, with each neuron now connected on each side to its r nearest neighbors, improved the synchronization process as assessed from a rise in the CI from 0.11 to >0.3 (Fig. 4A). Networks with larger connection radii produced oscillations of greater amplitude and lower frequency (Fig. 4B). The period at resonance was equal to approximately four times the mean delay d of the synaptic response (Fig. 3B).
Increasing the connection radius made the oscillations more robust at nonoptimal values of τ. As each neuron was connected to neighbors positioned at various distances, within a fixed axon radius r and with delays proportional to distance, a rectangular distribution of delays was formed. Inhibition through synapses with distributed delays made resonance less dependent on the appropriate value of τ because synchronously firing neurons at increasing distances activated their synapses with staggered delays. More particularly, in a linear array, the most distant afferent neuron would activate its synapse with a delay almost twice the mean delay, i.e., halfway in the oscillation cycle. This temporal summation of IPSCs induced a more lasting inhibition in the postsynaptic neuron, comparable with that induced by the activation of a single slowly decaying synapse. As a consequence, with a realistic decay time constant as small as 1.5-3 msec (Bartos et al., 2001; Carter and Regehr, 2002), resonance could be induced in the entire gamma frequency range and beyond by varying only the size of the axonal connection radius r (Figs. 3B, 4B).
The same relationship between mean delay and resonance was observed in two-dimensional networks, provided the connection probability between two neurons, or the connection weight, tapered off with distance to obtain a uniform distribution of synaptic delays or weighted delays. For example, the tuning curve in Figure 5A peaked at 78 Hz, i.e., at a period four times the mean weighted delay of network connections (3.2 msec). The tuning was sharp, oscillations being produced only within a narrow frequency band around fR for levels of excitation ranging from 12 to 375 spikes/sec, i.e., for all but the two rightmost excitation levels plotted in Figure 5B. Hence, although resonance was achieved when the neuronal firing rate equaled the oscillation frequency, oscillations at frequencies close to fR were generated over a broad domain of firing rates, encompassing almost the entire dynamic range of the neurons. Indeed, levels of excitation that produced average firing rates from 22 to 913 spikes/sec in a disconnected network evoked oscillations in the connected network restricted to a 47-149 Hz band (Fig. 5B).
Figure 5, C and D, illustrates in greater detail some observations already mentioned for the one-dimensional network. At resonance (Fig. 5C), the mean neuronal firing rate (73.4 spikes/sec) was close to the network frequency on the periodogram (78 Hz; Fig. 5C, a). The neurons fired a single spike at each cycle, producing a single peak at 13 msec on the averaged interspike interval histogram (Fig. 5C, b). The raster plot (Fig. 5C, c) and action potential trajectories (Fig. 5C, d) further demonstrate that synchrony was very robust though not very precise (CI = 0.287). At the next lower level of excitation (Fig. 5D), network frequency (72 Hz; Fig. 5D, a) was higher than the average firing rate (45 spikes/sec) and close to fR (83 Hz). The interspike interval histogram (Fig. 5D,b) showed a main peak at 25 msec, i.e., approximately twice the period at resonance. The synchrony induced at the onset of excitation rapidly faded away (Fig. 5D, c), but 72 Hz oscillations kept waxing and waning throughout the simulation (Fig. 5D, d).
Resonance arises in networks with sparse, asymmetrical connections
To determine the critical synaptic number for resonance, the mean number of connections per neuron was varied, whereas total synaptic weight was kept constant, in a network without topographic ordering, i.e., with each neuron able to connect to each other neuron with a single delay d = 1 msec. A network with a mean number of five synapses per neuron produced at resonance a power close to that achieved in the fully connected network (Fig. 6A, open circles). This threshold did not depend on the connection rule because power hardly increased when the network was randomly connected using the restriction that all neurons received the same number of synapses (Fig. 6A, filled diamonds). This absolute threshold for synchronization was further independent of the size of the network and the summed synaptic weight, and it was always less than the number of synapses needed to obtain synchrony in a network without delays (Fig. 6A, dashed curve).
The tuning curves (Fig. 6B) of the sparsely connected network (5 synapses per neuron; open squares; fR = 225 Hz) and the fully connected network (100 synapses; filled squares; fR = 214 Hz) were centered at approximately the same resonance frequency as the tuning curve of the network with only nearest neighbor coupling, shown in Figure 1F (fR = 210 Hz). The sparsely connected network was the more narrowly tuned (Fig. 6B) because antiphase synchronization produced a broader falloff at low frequencies in the fully connected network.
Effects of noise and heterogeneity
The randomization of synaptic weights and resting membrane potentials and the randomness in the activation of excitatory synapses made the present network heterogeneous and noisy.
We excited the nontopographic network in a noiseless manner by injecting a constant current into each neuron. The raster plots in Figure 6, C and D, compare the firing patterns of the sparsely (Fig. 6C) and fully connected (Fig. 6D) network during excitation by randomly firing fibers (Fig 6C, a, D, a) or an equivalent constant current (Fig. 6C, b, D, b). Removing the noise improved synchrony, as assessed from the increase in CI from 0.178 (Fig. 6C, a) to 0.381 (Fig 6C, b) and from 0.237 (Fig 6D, a) to 0.465 (Fig. 6D, b).
The residual asynchrony of the fully connected network (Fig. 6D, b) must be attributed to the heterogeneity of neurons and synapses. The homogeneous network produced a CI of 0.98 (data not shown). We quantified in an indirect way the present degree of heterogeneity by increasing the variability of current intensity among the neurons of the homogeneous network. Distributed currents drawn from a Gaussian distribution with a CV equal to 0.02 reduced the CI from 0.98 to 0.453, a value comparable with the CI of 0.465 in the raster plot of Figure 6D, b.
The noise level can also be changed by varying the number of afferent fibers. Because the membrane potential has a lower variance when large numbers of afferents excite a neuron through weak synapses, disconnected neurons fire more regularly when the number of fibers increases, whereas smaller numbers of fibers, activating stronger synapses, synchronize the disconnected neurons more precisely (Maex and De Schutter, 1998a,b; Maex et al., 2000). We examined which effect predominates in the generation of network oscillations and always observed that power increased with the size of the fiber population (Fig. 7A). In the limit of an infinite number of fibers, excitation by fibers is identical to the injection of a constant current. The tuning curve constructed by current injection did not differ in shape from that obtained with afferent fibers (Fig. 7B, filled squares; fR = 209 Hz).
Varying the degree of heterogeneity altered the peak power over almost two orders of magnitude without affecting much fR, although there was a tendency for low-frequency oscillations to be more robust than high-frequency oscillations.
Adding gap junctions stabilizes the oscillations but does not change fR
Closely spaced inhibitory neurons form dendritic gap junctions [>60% of neuron pairs separated <50 μm (Galarreta and Hestrin, 1999; Gibson et al., 1999) and all reciprocally connected pairs (Galarreta and Hestrin, 2002)]. In accordance with previous modeling studies (Traub et al., 2001; Bartos et al., 2002), the power of the network oscillations increased in a non-frequency-selective manner when adjacent neurons were connected by purely resistive electrical synapses (Fig. 7B, open diamonds). Hence, although the electrotonically transmitted afferent spike conceals the initial phase of the IPSP (see Tamás et al., 2000), this apparent increase in delay of the IPSP onset does not decrease the resonance frequency of the network (fR = 210 Hz).
Mechanism of synchronization and frequency tuning
Pairs of reciprocally connected inhibitory neurons tend to fire synchronously if each neuron, on firing, resets in an appropriate way the firing cycle of the paired neuron. In the present model neuron, activation of an inhibitory synapse delayed the generation of a spike, and the resulting increase of the interspike interval was greater the later in the firing cycle the synapse was activated (Fig. 8A, phase-response curve in B). As a consequence, of a pair of neurons reciprocally connected through inhibitory synapses, the neuron firing earlier received inhibition later and postponed its next spike more (Fig. 8C). Repetition of this mechanism progressively leads to a more precise synchronization of successive pairs of spikes (Ernst et al., 1995).
More precisely, let two mutually inhibitory neurons, denoted by superscripts 1 and 2, fire regularly with the same period T in a disconnected network. If the neurons generated their most recent spikes, indicated by subscript i, at ti1 and ti2 = ti1 + Δi, with Δi < d, then they receive inhibition in a connected network at time: 3 i.e., with a phase equal to: 4 If the phase-response curve has a linear form: then the neurons will generate their next spikes at time: 5 giving: 6 The neurons synchronize asymptotically if|Δi + 1| <|Δi|, which requires a <1 (assuming a > 0). Figure 8C depicts the sequence of events for phase-response curves with slope a <0.5. For 0.5 < a < 1.0, the paired neurons still converge to synchrony, but the order of their spikes reverses at each cycle.
The frequency tuning of the network oscillations and, more particularly, the 1/(4d) frequency observed for resonance can be understood from the effect of inhibition on the instantaneous firing probability of a neuron. When the first neuron of a pair has fired, the second neuron lowers its firing rate after the onset of inhibition, i.e., after an interval equal to the delay d. Alternatively, the second neuron may have fired first, preceding the first neuron within an interval of the same duration, so that the second neuron will preferentially fire in an interval [-d, d] around each spike of the first neuron (Fig. 8D). The sine wave best covering this 2d window of increased firing has a 4d period, and resonance occurs when both neurons, driven by the appropriate amount of excitatory input, fire on average with 4d intervals. In addition, for 4d to be the optimal period, inhibition must ensure that the 2d windows of increased firing alternate with 2d intervals of suppression of firing. This effect of the duration of inhibition explains the existence of an optimal decay time constant for each delay in Figure 3A.
Delayed reciprocal inhibition is able to synchronize, with zero phase lag, homogeneous networks of pulse-coupled oscillators and was proposed to be a mechanism of neuronal synchronization (van Vreeswijk et al., 1994; Ernst et al., 1995). Simulating heterogeneous hippocampal interneuron networks, Bartos et al. (2002) concluded that “a rapid inhibitory signal generated with a certain delay is a very effective synchronization signal.” The present study demonstrates that the resulting oscillations in network activity are limited to a narrow frequency band constrained by the IPSC latency but rather insensitive to the IPSC strength and decay time constant. Small, realistic axonal and synaptic delays, considered negligible compared with the low-pass time constant of the neuronal membrane (Manor et al., 1991), have a dramatic effect on the frequency spectrum (Fig. 1). The induced synchrony may be less precise than the synchrony that was achieved without explicit delays in more homogeneous networks (Wang and Buzsáki, 1996), but the robustness of the oscillations, their sharp tuning and frequency constancy, and their emergence in sparsely connected networks favor a model of delay-induced synchrony for fast brain rhythms. The extreme frequency control in these networks can be appreciated from the fact that variations in network heterogeneity did not affect the resonance frequency even when power changed by almost two orders of magnitude.
Oscillations obeying a 1/(4d) rule are observed in many biological systems (May, 1976; Glass and Mackey, 1988; MacDonald, 1989). In previous analytical models of oscillatory behavior in populations of neurons, delays attributable to signal propagation along axons and dendrites were lumped into the time constants of linear or nonlinear differential equations (Wilson and Cowan, 1972; Freeman, 1975). However, the incorporation of discrete delay variables can profoundly alter the dynamics of a system (MacDonald, 1989). Linear first-order delay systems can oscillate with a 4d period, and nonlinear delay differential equations can exhibit, at a critical delay value, a bifurcation from a steady state to a limit cycle solution with a period approximately equal to four times the delay period. In Appendix, we derive a formal explanation for our present findings.
Sources of delays in neural networks
The present results hold irrespective of the source of the delay between the timing of the action potential of an inhibitory neuron and the onset of the current response in the paired postsynaptic neuron. Because the range of the delays that neural circuits are able to produce constrains the frequency range of the oscillations to which the present model may apply, we briefly discuss some prevalent mechanisms of delayed inhibition.
Although inhibitory neurons are fast-processing units, with an axon tree often confined within 600 μm from the soma (Buhl et al., 1994; Sik et al., 1995), some (partly myelinated) axons can extend several millimeters in hippocampus (Sik et al., 1994) and neocortex (Kisvarday et al., 1993; Gupta et al., 2000). Salin and Prince (1996) electrically evoked monosynaptic IPSCs with latencies of 1 to >6 msec in neocortical pyramidal cells in vitro and derived slow speeds for spike propagation (0.06-0.2 m/sec; mean, 0.1 m/sec). From our simulation data, a mean latency of 4 msec induces resonance at 62.5 Hz (Fig. 3), which is within the frequency range of the oscillations recorded in visual (Gray and Viana Di Prisco, 1997) and auditory neocortex (Brosch et al., 2002).
The reported latencies of unitary IPSCs, evoked by firing an impaled presynaptic neuron, are usually much shorter (mean, 0.8 msec for pairs <200 μm apart in dentate gyrus; Bartos et al., 2001), but paired recordings may be biased to small interelectrode distances at which the probability of finding connected neurons is highest. IPSP latencies of 3.2-8.6 msec (mean, 5.4 msec) were measured for pairs 153-445 μm apart in striatum (Tunstall et al., 2002), and in hippocampus, unitary IPSCs between lacunosum-moleculare interneurons and pyramidal cells showed latencies, albeit at room temperature, of 2.4-7.2 msec (mean, 4.2 msec; Bertrand and Lacaille, 2001).
In addition, it is the mean distance to the entire set of postsynaptic neurons that determines the resonance frequency. This mean distance can be estimated from the experimentally derived connection probability function, provided the probability at each distance is corrected for the varying numbers of neurons available. For axons with a disk- or sphere-like arborization, the number of candidate target neurons increases linearly or quadratically with distance from the source neuron. Hence, if the connection probability is observed to fall off according to an exponential function with space constant δ, the actual mean distance increases from δ to 2δ and 3δ for axons branching in two and three dimensions, respectively.
The present mechanism is not restricted to purely inhibitory networks. The 1/(4d) rule generalizes to circuits with intercalated excitatory neurons (for example I1 → E1 → I2 → E2 → I1 instead of I1 → I2 → I1). Here resonance is predicted to arise at an oscillation period equal to four times the delay of the disynaptic response I1 → E1 → I2. If the response delay is much greater for excitatory than for inhibitory connections, then the circuit will resonate with a period approximately four times the delay of excitation, and, consequently, the excitatory neurons will lead the inhibitory neurons by one-fourth of a cycle. Such a phase relationship has been observed in olfactory systems (Freeman, 1975; Bazhenov et al., 2001) and hippocampus (Csicsvari et al., 2003).
Finally, circuits with excitatory feedback can generate trains of spikes with delays of tens to hundreds of milliseconds. In visual cortex, feedback excitation was proposed to provide the long delays of inhibition needed to generate directionally selective responses to slowly moving stimuli (down to <1 Hz; Maex and Orban, 1996). It is noteworthy in this respect that the property of directional selectivity, which can be implemented in a circuit with lateral inhibition, exhibits the same f = 1/(4d) relationship, f being the temporal frequency of the moving stimulus evoking maximally selective responses (van Santen and Sperling, 1985).
The present results and the above arguments lead to the following predictions, some of which are substantiated by recent findings. The delay of monosynaptic or polysynaptic inhibition is the critical parameter determining the resonance frequency of an inhibitory network. Hence, oscillations with different characteristic frequencies are expected to be generated by microcircuits involving different types of interneurons (Klausberger et al., 2003). In these interneurons and in their target projection neurons, IPSCs are predicted to be evoked with variable latencies, centered at approximately one-fourth of the oscillation period (Traub et al., 1996). Because the multitude of microcircuits in the brain could generate a continuum of delays, probably with slight differences emerging between different areas or developmental stages, oscillations in the nervous system could form a continuum rather than being divided into a few discrete frequency bands (Csicsvari et al., 1999; Grenier et al., 2001). Interneurons involved in circuits with low characteristic frequencies are predicted to fire multispike bursts during each cycle. Finally, although the delay of reciprocal inhibition might not be very amenable to experimental manipulation, changing the size of the stimulus would vary the effective mean delay. Stimuli exciting a focus smaller than the connection radius of an inhibitory network do not recruit long-distance neuron pairs, leading to faster oscillations than predicted from the mean delay of the circuit (Bartos et al., 2002, their Fig. 4E). Fast synchronous oscillations are therefore predicted to be more narrowly localized than low-frequency oscillations, although long-distance excitatory connections may contribute to this difference (Ermentrout and Kopell, 1998; Pauluis et al., 1999; Kopell et al., 2000).
We propose that resonant synchronization induced by the delay of (monosynaptic or polysynaptic) inhibitory connections contributes to the emergence and frequency tuning of all types of oscillations in which inhibitory synapses are involved, and that delay-induced synchronization should be considered for fast oscillations in particular (>40 Hz). Some types of fast oscillations appear to persist, however, in the absence of synaptic transmission (200 Hz “ripples”), and axoaxonal gap junctions may be essential for their generation (Schmitz et al., 2001). For other types of oscillations, additional tuning mechanisms may be important, such as intrinsic neuronal resonance and synaptic dynamics (Gupta et al., 2000; Beierlein et al., 2003), which were not included in the present model.
Let A(t) be the population activity at time t, F the neuronal transfer function, I the constant level of external input, K(t) the synaptic impulse response, and d the combined axonal and synaptic delay. Then the dynamics of a fully connected, homogeneous network can be described by: A1 where · denotes convolution in time (Gerstner, 2000). During asynchronous firing, the population activity is considered constant over time: A(t) = C = F(I - C). Linearization about the asynchronous state yields: A2 The linearized equation, after substituting y(t) for A(t) - F(I - C) and a for the derivative of F is: A3 We take K(t) = (1/τ) exp(-t/τ) and follow the reasoning developed by MacDonald (1989), looking for a periodic solution y = exp(i ω t), with ω = 2π/T. We consider two cases. The first case explains the observation in Figure 3B that for networks with low resonance frequency, fR is slightly >1/(4d) but almost completely independent of the decay time constant τ. The second case matches the observation that networks with a high resonance frequency can have an fR of ≤1/(4d), which in addition is more dependent on the value of τ.
At low resonance frequencies (case 1), the synaptic kinetics is much faster than both the membrane time constant λ and the oscillation period T, so that the delayed feedback can be considered to be approximately instantaneous. This reduces Equation A3 to: A4 The resulting transcendental characteristic equation: A5 has real roots for cos(ω d) = -1/a and sin(ω d) = ω λ/a, or: A6 Provided a > 1, solving for ω d yields π/2 < ω d < π, giving 2 d < T <4 d. Hence, fR > 1/(4d) but independent of τ.
At high resonance frequencies and hence high levels of excitation (case 2), the membrane potential of a neuron is strongly depolarized, and voltage deflections are small (although sufficient to pass the firing threshold) because of the low-pass characteristics of the membrane. Hence the contribution of the leak current to membrane dynamics can, approximately, be ignored. Equation A3 now reduces to: A7 The resulting transcendental characteristic equation is: A8 For very fast synapses in the limit τ → 0, the solution is ω d = π/2, or T = 4d. For τ > 0, the present case is an instance of a distributed delay (the synaptic kernel) with a gap (the synaptic delay d), and T >4d with T increasing for increasing τ (MacDonald, 1989).
Note added in proof. It was brought to our attention that Brunel and Wang (2003) reported similar results recently.
This work was supported by Inter-University Attraction Pole Programme Grant 5104, Fonds voor Wetenschappelijk Onderzoek Grant G.0401.00 (Vlaanderen), and European Union Grant QLRT-2001-01241. We thank Mike Wijnants for graphics expertise and anonymous reviewers for contributions, particularly to Appendix.
Correspondence should be addressed to Reinoud Maex, Born-Bunge Foundation, University of Antwerp, Universiteitsplein 1, B-2610 Antwerp, Belgium. E-mail:.
Copyright © 2003 Society for Neuroscience 0270-6474/03/2310503-12$15.00/0