Abstract
Slow population activities (SPAs) exist in the brain and have frequencies below ∼5 Hz. Despite SPAs being prominent in several cortical areas and serving many putative functions, their mechanisms are not well understood. We studied a specific type of in vitro GABAergic, inhibition-based SPA exhibited by C57BL/6 murine hippocampus. We used a multipronged approach consisting of experiment, simulation, and mathematical analyses to uncover mechanisms responsible for hippocampal SPAs. Our results show that hippocampal SPAs are an emergent phenomenon in which the “slowness” of the network is due to interactions between synaptic and cellular characteristics of individual fast-spiking, inhibitory interneurons. Our simulations quantify characteristics underlying hippocampal SPAs. In particular, for hippocampal SPAs to occur, we predict that individual fast-spiking interneurons should have frequency–current (f–I) curves that exhibit a suitably sized kink where the slope of the curve decreases more abruptly in the gamma frequency range with increasing current. We also predict that these interneurons should be well connected with one another. Our mathematical analyses show that the combination of synaptic and intrinsic conditions, as predicted by our simulations, promotes network multistability. Population slow timescales occur when excitatory fluctuations drive the network between different stable network firing states. Since many of the parameters we use are extracted from experiments and subsequent measurements of experimental f–I curves of fast-spiking interneurons exhibit characteristics as predicted, we propose that our network models capture a fundamental operating mechanism in biological hippocampal networks.
Introduction
Despite constituting only 10–15% of the total neuronal population, the inhibitory population of interneurons forms a critical component of the mammalian hippocampus underlying a wide range of population activities. Notably, it generates gamma rhythms (25–100 Hz) that are thought to be responsible for spatial information and representation of memory (Buzsáki, 2006). Experimental and modeling studies (Wang and Buzsáki, 1996; Vida et al., 2006; Bartos et al., 2007; Sohal et al., 2009) ascribe the ability of the inhibitory population to support gamma rhythms to the highly connected network of fast-spiking interneurons via fast GABAergic synapses. The fast kinetics of GABAergic synapses and the high interconnectivity allow synchronization among fast-spiking interneurons, thus creating high-frequency population activities. However, inhibitory networks also sustain slow population activities (SPAs) of much lower frequencies (<5 Hz) (Schwartzkroin and Haglund, 1986; Straub et al., 2000; Papatheodoropoulos and Kostopoulos, 2002; Wu et al., 2005b).
In particular, Wu et al. (2002, 2005b) have established that coherent activity of inhibitory interneurons (including fast-spiking interneurons) underlies SPAs in mouse hippocampus. Given the disparity between the timescales of individual cell firings (∼100 Hz) and SPA network frequencies (<5 Hz), how do we reconcile the apparent paradox that these fast-spiking interneurons with fast synaptic kinetics can support gamma rhythms as well as SPAs? This timescale paradox is also inherent in many brain population activities (Bullmore et al., 2009; He et al., 2010). Thus, its resolution provides insight into the complexity of brain dynamics, where rhythms of frequencies spanning four orders of magnitude routinely coexist (Buzsáki and Draguhn, 2004).
In this article, we focus on the in vitro hippocampal SPAs (Fig. 1A) of Wu et al. (2005b) to explore how both intrinsic and network characteristics shape slow activities in networks of fast-spiking interneurons. A system can have emergent population slow timescales during a phase transition (where a network changes global firing patterns) (Landau and Lifshitz, 1980; Plischke and Bergersen, 2006). We have previously extracted synaptic constraints during the ongoing population activity (Ho et al., 2009). Based on these constraints, we developed interneuronal network models and found that SPAs can be generated in model networks of fast-spiking interneurons with rapid synaptic properties. From our simulations, we uncovered the relationship between intrinsic cellular characteristics of constituent interneurons and the emergence of SPAs in inhibitory networks. Specifically, individual interneurons must have frequency–current (f–I) curves with an abrupt slope decrease as I increases, a “kink,” for the inhibitory network to sustain SPAs. Using mathematical analyses, we showed that such f–I curves promoted the coexistence of multiple network firing states (i.e., network multistability). Population slow timescales are the result of the network switching from one state to another due to excitatory synaptic fluctuations. Furthermore, we experimentally determined that these f–I characteristics, as required for the existence of SPAs, are present in CA3 fast-spiking cells where SPAs robustly occur.
SPAs are also prominent in cortical circuitries (Steriade et al., 1993; Jarosiewicz et al., 2002) other than hippocampal networks. Cortical SPAs are thought to serve functions including memory consolidation (Buzsáki, 1989) and synaptic plasticity (Tsukamoto-Yasui et al., 2007). Their mechanisms are, however, not well understood in terms of excitatory and inhibitory interactions and balances. Although our work is on GABAergic, inhibition-based SPAs in hippocampus, our findings may generalize to other circuitries. Our results shed new light on the synergistic effects between intrinsic and synaptic properties underlying SPAs.
Materials and Methods
Experiments
Field potential recordings in hippocampal slices.
Hippocampal slices (thickness, 1 mm) were prepared from 1- to 3-month-old C57BL/6 mice of either sex, as previously described (Wu et al., 2005a,b). All animal procedures were performed in agreement with the policies and guidelines of the Canadian Council on Animal Care. To achieve sufficient in vitro perfusion and oxygenation of the 1 mm thick slices, we separated the dentate gyrus from the CA1 area while keeping functional connections of the dentate gyrus-CA3 projection. The slices were perfused with an artificial CSF (ACSF) at 32–33°C. The ACSF contained the following (in mm): 3.5 KCl, 1.25 NaH2PO4, 125 NaCl, 25 NaHCO3, 2 CaCl2, 1.3 MgSO4, and 10 glucose. Extracellular recordings were made with glass pipettes filled with a solution containing 200 mm NaCl and 5 mm HEPES, pH adjusted to 7.4. A CA3 extracellular recording electrode was placed largely in the somatic or cell body layer.
Spontaneous SPAs were observed to arise in the CA3 area (Fig. 1A) and to spread to subicular and entorhinal cortical areas (Wu et al., 2006), and simultaneous cellular recordings indicated that coherent activities in inhibitory cells underlay these SPAs (Wu et al., 2002, 2005b). These hippocampal SPAs were also previously referred to as spontaneous rhythmic field potentials in Wu et al. (2005b) and as basal sharp waves in Ho et al. (2009).
Whole-cell recordings of GABAergic inhibitory interneurons.
Acute hippocampal slices (thickness, 300 μm) were cut from 18- to-21-day-old Wistar rats of either sex using a VT1200 vibratome (Leica), as described previously (Sauer and Bartos, 2010). All animal procedures were performed in agreement with national and institutional guidelines. For storage and superfusion of slices, a physiological solution was used containing the following (in mm): 125 NaCl, 25 NaHCO3, 2.5 KCl, 1.25 NaH2PO4, 25 glucose, 2 CaCl2, 1 MgCl2. Slices were incubated at 34°C for 15–30 min and subsequently stored at room temperature. Solutions were equilibrated with a 95% O2–5% CO2 gas mixture throughout the procedure. Patch pipettes for whole-cell patch-clamp recordings were pulled from borosilicate glass capillaries (outer diameter, 2.0 mm; inner diameter, 1.0 mm) and had a resistance of 2.5 to 4.0 MΩ. Pipettes were filled with K-gluconate-based internal solution as follows (in mm): 120 K-gluconate, 20 KCl, 10 EGTA, 2 MgCl2, 2 Na2ATP, 10 HEPES).
Interneurons with somata located in or close to the principal cell layer of CA3 were recorded under visual control using infrared differential interference contrast videomicroscopy (Doischer et al., 2008). Whole-cell recordings were made using an Axopatch 200B amplifier (Molecular Devices). The resting membrane potential was determined after breakthrough in the I = 0 mode and was −65 to −55 mV. Series resistance (Rs, 9–25 MΩ) was monitored and compensated throughout the experiment (100%, time lag 1 ms). Signals were filtered at 10–20 kHz (4-pole, low-pass Bessel filter) and digitized at 20–40 kHz using a Power 1401 digital interface (Cambridge Electronic Design). Pulse generation and data acquisition were performed using FPulse (U. Fröbe, Physiological Institute, University of Freiburg, Freiburg, Germany) running under Igor Pro 5.01 on a personal computer. For data analysis, Stimfit (C. Schmidt-Hieber, University College London, London, UK) and MatLab (Mathworks, http://www.mathworks.com) were used. All recordings were performed at 32–34°C.
To determine f–I curves (see Fig. 6B,C), cells were held in the current clamp mode and 1 s long depolarizing current injections are applied. To obtain a high f–I resolution, amplitudes were increased incrementally with step sizes of 5–10 pA. The maximal current injection was set as the amplitude reaching maximal discharge frequency. Interneurons were classified as fast-spiking (FS) on the basis of maximal discharge frequency >100 Hz, and an adaptation ratio >0.5 was determined from the ratio of the first and the last interspike intervals during a train of action potentials at ∼100 Hz (1 s).
We note that the recorded maximal discharge frequencies were comparable to those previously determined in FS interneurons in CA3 (Bartos et al., 2002) and dentate gyrus (DG) (Doischer et al., 2008) of mice. This indicates that f–I relations are independent of the particular hippocampal area and of whether rats or mice are used.
Detection and quantification of kinks in f–I curves.
We analyzed f–I curves of FS interneurons to determine whether they exhibited a kink, defined as an abrupt decrease in slope. To determine the position of the most prominent kink in the recorded f–I curves, the local slope (i.e., Δf/ΔI) of the f–I curve was estimated by linear fitting of short sections (length, 30 pA, analysis resolution) of the f–I curve in question. Figure 6, B and C, shows examples of
To directly compare kinks in f–I curves from experiments with those in our models, we used the same protocol described to quantify the kink characteristics of our model interneuron f–I curves (see Fig. 3A,B, analysis resolution: 0.015 μA/cm2). However, to apply the same protocol, we had to approximate the experimental situation in our model f–I curves. We did this by including a very small amount of noise (i.e., ge0 = 0 mS/cm2, σe = 0.00027 mS/cm2; see equation set 4 below) to our model interneurons for kink calculations. In Table 1 we show the results of our analyses of six model interneuron f–I curves.
Summary of model interneuron f–I curves and mean-field results
Simulations
Model description and parameter choices.
Each inhibitory interneuron was modeled by the following set of equations (1 ≤ i ≤ N = 120 denotes cell number, V is membrane potential in millivolts, and t is time in milliseconds),
where Iinti represents the intrinsic current of the interneuron and is adapted from a simplified model by Izhikevich (2003, 2007). The specific capacitance C is assumed to be 1μF/cm2. Iinti has the following form,
Specifically, Vth and Ith represent, respectively, the values of membrane potential and external current when the neuron is at threshold, α controls the increase in the firing rate of the interneuron just above threshold, while p controls the functional form of this increase expressed by the model interneuron. For example, p = 2 represents a model interneuron with a saddle node on an invariant circle type bifurcation at the onset of firing, and p = 1 represents an integrate-and-fire type of model interneuron. We note that equation set 1 does not describe the suprathreshold dynamics for the case p = 1. As a result, for p = 1 we artificially made the interneuron fire and set Vi > Vmax whenever Vi crossed Vth from below. The variable ni is an auxiliary current for controlling the refractory period of the interneuron; it primarily affects the interneuron's firing behavior far from the threshold. Equation set 2 describes the behavior of this current. It resets to nreset whenever the neuron fires and decays exponentially with constant λ between individual spikes. Iext is the external current for the model interneurons. In single-cell cases (i.e., when all the model interneurons were disconnected from each other), Iext represents the injected current for the model interneurons. In network cases, the variable represents the “virtual” excitation of the model interneurons received from the excitatory cell population (equation set 4).
We used the above mathematical form for intrinsic current (equation set 2) and not the more traditional conductance-based model because there is a direct correspondence between the parameters of equation set 2 and the input–output characteristics of the neuron around threshold. Our models are capable of firing from 0 to up to ∼150 Hz so that characteristics of experimental fast-spiking interneurons could be adequately represented. Figure 3 shows sample f–I curves and sample output of voltage versus time from a simulation of a single firing model interneuron under constant injected current.
The synaptic current Isyni is represented as,
where Esyn is the inhibitory reversal potential and the si (i = 1–120) represents the synaptic gating variables of the neurons that are synaptically connected (with strength gsyn) to the neuron in question. For computational efficiency, we used a discontinuous model to represent the opening of the synapses (equation set 3) whenever there was a spike (comparative simulations indicate that results are similar to those using a continuous model). Parameters τsyn and smax denote, respectively, the inhibitory synaptic decay time constant and the value of the maximal increased opening of the synaptic gates per spike.
To mimic the excitation that the inhibitory population receives from the pyramidal cells, Iext has the form (Rudolph et al., 2004)
where Ee is the excitatory reversal potential and τe is the excitatory time constant. The excitatory conductance ge is a stochastic variable following the Ornstein-Uhlenbeck process (equation set 4) with mean ge0 and SD σe (excitatory fluctuations). χe underlies the stochasticity of ge. To first order, ∫t1t2 χe(s)ds, t2 − t1 > 0 can be, but does not necessarily have to be, represented as N(0,t2 − t1), where N(0,t2 − t1) is Gaussian distributed random variable of mean 0 and variance t2 − t1 (Greiner et al., 1987; Gardiner, 2004; Seydel, 2004).
Therefore, our computational network model consisted of explicit mathematical representations for inhibitory cells (equation sets 1–3) and virtual representations for the excitatory cell population (equation set 4). Each model inhibitory interneuron was connected with all the other inhibitory interneurons in the network (unless otherwise stated) while receiving noisy excitation from the virtual excitatory cell population. We changed the intrinsic characteristics (f–I curves) of inhibitory interneurons by adjusting the values of α, p, λ, and nreset.
We used moderately sized networks of 120 fast-spiking cells. We consider this size to be large enough as it has been successfully used by others to obtain mechanistic insights (for review, see Wang, 2010). Given that the hippocampal SPAs are due to coherent inhibitory cell activities and do not seem to preferentially involve specific inhibitory cell types (Wu et al., 2002, 2005b), we reasonably focused on fast-spiking cells in our network models since they constituted the majority of inhibitory cells (Vida, 2010). We used an average of the summed synaptic variables (i.e., Σi si/N) to represent the local field potentials (LFPs). This is reasonable given that LFP generation is mainly due to synaptic currents. Earlier simulations using synaptic currents rather than synaptic variables showed that this did not affect our conclusions. Since it was easier to use synaptic variables, we did so.
We divided our simulations into different sets. Each set of simulations consisted of 828 independent simulations of the 120-cell inhibitory network for 65 s (for homogeneous simulations only; for nonhomogeneous simulations, the simulated time was 8–10 s) with different gsyn and σe values (gsyn values from 0.2 to 7.0 × 10−2 mS/cm2 with increment 0.1 × 10−2 mS/cm2, and σe values from 0.027 to 0.302 × 10−2 mS/cm2 with increment 0.025 × 10−2 mS/cm2). The first 5 s of any homogeneous simulation were excluded from further analysis. The different simulations within each set of simulations were characterized by having identical intrinsic properties (i.e., f–I curves) of inhibitory interneurons but different synaptic properties between the interneurons. In other words, between any two simulations within the same set, the intrinsic properties of, say, model interneuron 1 of one simulation were the same as those of model interneuron 1 of the other simulation, and so on for model interneurons 2–120. However, the synaptic properties (e.g., gsyn) between, say, model interneurons 1 and 2 in one simulation were generally different from those of the other simulation. The gsyn values were drawn from data of unitary inhibitory conductance measurements from Bartos et al. (2002). The σe values were chosen to incorporate the range of values estimated directly from experimental data of our in vitro hippocampal SPAs experiments (Ho et al., 2009). Table 2 shows the common parameter values we used across all sets of simulations, unless specified otherwise. Table 3 shows other simulation parameters that were changed across different simulations to vary the synaptic and intrinsic conditions of the simulated networks.
Constant parameters used for the inhibitory network simulations
Other parameters used for the inhibitory network simulations
To ensure that our simulated networks received similar excitation when making comparisons, we made the constructed f–I curves have the same threshold and coincide with each other at exactly one point above the threshold. Figure 3 shows examples of f–I curves of model inhibitory interneurons we used for simulation. The ge0 (mean excitation received by model inhibitory interneurons; see equation set 4) value was common across all sets of simulations unless otherwise specified. The value of ge0 was chosen so that neurons from different sets of simulations would fire at approximately the same frequency were they synaptically isolated from inhibitory interactions. We note that this value is in line with estimates from Ho et al. (2009).
We designed custom C++/Perl codes for the simulations. Owing to the enormous amount of random numbers for the stochastic integration (χe of equation set 4), a robust random number generator was required. We used a pseudorandom number generation algorithm (ran2), as described by Press et al. (1997), which is known to be robust and to have a very long period (>2 × 1018). For a network of 120 neurons with integration time step of 0.01 ms, this length of period enables simulation up to at least 1013 ms (1013 × 120 × 100 = 1.2 × 1017, which is at most 6% of the period of the random number generator). We used the Box–Muller algorithm (Press et al., 1997) to transform uniformly distributed random variables (i.e., outputs of ran2) into Gaussian distributed random variables. We performed most of the simulations on the GPC supercomputer at the SciNet HPC Consortium (Loken et al., 2010) (http://www.scinethpc.ca/). We used the Euler algorithm for numerical integration. Integration time step for each simulation ranges from 0.001 to 0.01 ms. Preliminary simulations using a smaller number of model interneurons were performed using XPPAUT (Ermentrout, 2002) (http://www.math.pitt.edu/∼bard/xpp/xpp.html).
Metric and frequency calculation for simulated SPAs.
To determine the occurrence of SPAs and to quantify the strength and frequencies of SPAs, we devised a metric based on the power spectra (i.e., discrete Fourier transform) of the simulated network outputs (Σ si/N). The main idea was the recognition of a “signature” pattern in the power spectra whenever the outputs exhibited slow population activities. We calculated the spectra using the entire dataset after discarding the first 5 s of transients. Spectra were normalized by the original signal and the total time. We used the subroutines in the FFTW (Fastest Fourier Transform in the West) package (Frigo and Johnson, 2005) (http://www.fftw.org) to compute discrete Fourier transform for the power spectra and metric.
Spectra examples can be seen in Figures 1B and 2. One can immediately notice that there are two “humps” in each of the spectra from 0 to 100 Hz, one toward the low end of the frequency and the other peaking at ∼40 Hz. The peak toward the low frequency (<5 Hz) was responsible for the slow timescale of the simulated SPAs while the peak at ∼40 Hz was due to the individual firing of inhibitory interneurons. We visually inspected all simulations and their corresponding spectra and noted that the characteristic “two-hump” shape of the power spectrum was a good indicator for the presence of SPAs in simulations. We thus developed a heuristic approach to quantify simulated SPAs based on the shape of the power spectrum. For any simulation, we regarded SPAs as present whenever (1) the peak between 0.1 Hz and 5 Hz (plow) was at least 2× the average power between 0.1 and 5 Hz, (2) the peak above 5 Hz (phigh) was at least 2× the average power density above 5 Hz and above some threshold (phigh,t), and (3) plow was at least 0.4× phigh. We used (phigh,t) = 0.05 × 10−5 Hz−1. We note that every power spectrum obtained from simulation was produced by discrete Fourier transform from the original simulation with a sampling interval of 0.02 ms, regardless of the size of integration time steps used in the original simulation. Thus, for the strength of simulated SPAs (a dimensionless quantity), we used the value of plow divided by the sampling interval once the above criteria (1, 2, and 3) were satisfied. These criteria deserve some justification. The first criterion was to ensure there was a respectable peak at the low end of the power spectrum, and the second criterion served the same purpose for the peak at the higher end of the power spectrum. The third criterion was to prevent the peak toward the low end of spectrum from being disproportionally small as compared to the peak at the higher end of spectrum. We computed the frequencies of simulated SPAs in the following way: for each simulation in which SPAs were defined to exist, we rearranged the power spectrum between 0.1 and 5 Hz in descending order of power density magnitude (which was possible because the spectra were discrete). To determine the frequencies of simulated SPAs, we included all the frequencies whose power density magnitudes were at least 80% of plow. We calculated mean and SD and plotted the results as error bars (see Fig. 7, bottom).
Thus, SPAs were deemed to have occurred in a particular network if the criteria for the spectra were satisfied, with the strength of SPAs calculated from the values of the low-frequency peaks (0.1–5 Hz) of the power spectra. The location of these low-frequency peaks in the spectra determined the dominant frequencies of the resulting network activities.
Mean-field equations.
To explain the observations from our network simulations, we approximated our network with a mean-field model. The dimensions of the mean-field model were much reduced from the full network model so that it would be amenable to mathematical analyses. Our mean-field model was developed from a rate-based formulation of neuronal dynamics by Amit and Tsodyks (1991) and Bressloff and Coombes (2000). In Appendix C of Ho (2011) (http://hdl.handle.net/1807/30056), a step-by-step reduction of the full simulation into the mean-field equations (equation set 5) has been outlined. The mean-field equations are an approximation of the full simulation in that they only capture the lowest order terms of the full simulation. In other words, fluctuation effects of the full simulation are neglected in the mean-field equations. We reproduce here the resulting equations for mathematical analyses,
where the parameters are the same as those in equation sets 1, 2, and 3. The above set of mean-field equations (i = 1 − 120) governs the temporally coarse-grained progression of the total synaptic current onto inhibitory cell i, given a homogeneous all-to-all network configuration with mutual inhibitory conductance of gsyn. f(…) is the f–I function of individual inhibitory interneurons. 〈V〉 is the average value of (Esyn − Vi) across all inhibitory interneurons and times. To further reduce the dimensionality of the set of 120 mean-field equations (equation set 5), we assumed the network was spontaneously broken into two clusters, a non-firing cluster a and a firing cluster b (simplest possible assumption). Once this assumption was made, the total number of mean-field equations was reduced to three,
where Na and Nb are the number of inhibitory interneurons in clusters a and b, respectively. Our aim was to find the ranges of values of gsyn, which yield stable solutions to equation set 6 for different sets of Na and Nb values. These sets of Na and Nb values are constrained by the total number of model inhibitory interneurons N. Since equation set 5 is temporally coarse grained, stable fixed points of equation set 6 (i.e., stable solutions for
Representing f–I characteristics analytically.
We represented the f–I characteristics of model inhibitory interneurons with analytical functions when performing mathematical analyses via mean-field equations. The interspike intervals were obtained analytically by first integrating equation set 1 with Isyni = 0 and nreset = 0. We call this first integrated result T1. To take into account the effects of the refractory periods caused by nonzero nreset, we added T1 to a constant time interval T2 to form the total interspike interval T. We note that, while we analytically derived T1 given the parameters α, p, Vmax, and Vreset, T2 was obtained by fitting
Results
SPAs can arise in model networks of fast-spiking inhibitory cells
To investigate the mechanisms leading to the emergence of SPAs in cortical circuits, we designed a network model composed of mutually connected fast-spiking inhibitory interneurons. In many of our network simulations implementing different interneuron single-cell models, inhibitory synaptic strengths, and fluctuation levels, we obtained SPAs characterized by periods of high network activity interspersed with low network activities at frequencies <5 Hz. Moreover, a subcollection of our simulated SPAs possessed characteristics similar to those of the in vitro hippocampal SPAs described by Wu et al. (2005b). These characteristics were as follows: first, the dominant SPA frequency occurred between 0.5 and 4.5 Hz; second, the duration of the high activity SPA episode was ∼200 ms; third, the firing frequency of model cells during SPAs was ∼40–80 Hz, as seen experimentally. Importantly, there was at least a twofold increase in the inhibitory conductance from low to high activity states in accordance with results from analyses obtained directly from the experimental data (Ho et al., 2009).
In Figure 1A, we show an extracellular field recording (LFP output) from a thick slice preparation exhibiting hippocampal SPAs. One representative network simulation exhibiting SPAs is shown in Figure 1B. We plot the average summed inhibitory synaptic variables (Σi si/N) as a representation for the local field potential (LFP analog) on the top of Figure 1B. In Figure 1, A and B, the bottom two panels are power spectra of the top panel (in log-log scale and regular scale, respectively). Although the LFP analog (Fig. 1B) is only a rough approximation to actual LFP (Fig. 1A), a clear similarity with the experimental data can be seen. Both simulation and experiment express slow activities that can be observed in the LFPs themselves on the top as well as their spectra on the bottom. An additional peak in the 40 Hz region due to the firing rate of individual fast-spiking cells is apparent in the model spectra (Fig. 1B).
SPAs in the experiment and the model have similar characteristics. A, CA3 experimental output from a thick hippocampal slice showing SPAs. The data have been bandpass filtered between 0.5 and 250 Hz. B, An example of simulated SPAs that resemble experimental data. Network output is shown as average summed inhibitory synaptic variables (∑ si/N) as an analog for LFP output. Parameters as in Table 2 with gsyn = 0.048, σe = 0.00152, α = 1, p = 2, λ = 0.81, nreset = 4. See Table 3 for units of parameters used. A, B, The power spectrum of the corresponding trace is shown on each lower-right panel, with the lower-left panel showing the spectrum on a log–log plot.
Cellular f–I characteristics control the occurrence of SPAs
We show two sample simulated SPAs in Figure 2. These two example outputs (Fig. 2A,B) were drawn from two different sets of simulations, which means that the model inhibitory interneurons used for the two simulations had different intrinsic properties (see Materials and Methods). The frequencies of these simulated SPAs can range from much less than 1 Hz to 5.0 Hz, as can be seen in the spectra in the bottom plots of Figure 2. From our simulations, we find that the increase in output synaptic conductances from the low to the high activity state can be as small as <2-fold to as large as 10-fold (Fig. 2A,B, compare the top plots). In a 15 ms time window, approximately one-tenth to one-third of the neurons in the network participate in firing during high activity states.
Model SPAs have a range of characteristics. A, B, Two examples of SPAs obtained in the network models. Parameters as in Table 2 with (A) gsyn = 0.040, σe = 0.00102, α = 0.41, p = 2, λ = 20, nreset = 4 and (B) gsyn = 0.044, σe = 0.00127, α = 1, p = 2, λ = 0.81, nreset = 4. See Table 3 for units of the parameters used above. A, B, Top, Temporal progression of the average synaptic output (Σsi/N) of the simulated network. We used (Σ si/N) as our LFP analog. Second from top, Raster plot of all the spiking neurons in the 120-cell network. Second from bottom, Intracellular trace of one of the neurons in the 120-cell network. Bottom, Power spectra of the LFP analog of the 65 s simulation. We excluded the first 5 s from analysis.
To quantify the strength of SPAs, we used a metric based on the power spectrum of the simulated SPAs (LFP representation). Briefly, using the model LFP power spectrum, a normalized version of the low-frequency height gives a value that quantifies the SPA strength from each simulation. This value is color coded in the plots of Figures 4, 5, and 7. The pseudo-color temperature indicates the strength of the SPAs. Two samples of strong SPAs are given by the green triangles in Figures 4A and 5A, and a sample of weak SPAs is given by the green diamond in Figure 4B. See Materials and Methods for further details of the metric computation.
f–I curves of model interneurons and sample output. A, f–I curves. Curves have Vth, Ith, Vreset, and Vmax parameters as in Table 2 with curve i, α = 0.41, p = 2, λ = 20, nreset = 4; curve ii, α = 1, p = 2, λ = 0.81, nreset = 4; curve iii, α = 10, p = 2, λ = 0.74, nreset = 32. B, f–I curves. Curves have Vth, Ith, Vreset, and Vmax parameters as in Table 2 with curve i, α = 0.2, p = 4, λ = 30, nreset = 4; curve ii, α = 1, p = 4, λ = 0.98, nreset = 4; curve iii, α = 3.5, p = 4, λ = 0.69, nreset = 4. C, Voltage versus time for a firing model interneuron. Parameters as in Table 2 with α = 1, p = 2, λ = 0.81, nreset = 4, Iext = 0.03 μA/cm2. See Table 3 for units of the parameters used above.
Dependence of SPAs on f–I curves for p = 2. Summary of three sets of simulations. A–C, Each set of simulations uses different intrinsic properties of model inhibitory interneurons while keeping p = 2. A, Set of simulations using α = 0.41, λ = 20, nreset = 4. B, Set of simulations using α = 1, λ = 0.81, nreset = 4. C, Set of simulations using α = 10, λ = 0.74, nreset = 32. See Table 3 for units of parameters used above. A–C, Top, Four selected simulations (LFP representations). For clarity of presentation, only a 30 s segment of the entire 65 s simulation is shown for each simulation. Middle, Summary of occurrence of simulated SPAs across synaptic parameter regimes (gsyn, σe). Warmer color (reddish-yellow) represents stronger SPAs. Black represents the absence of SPAs for the particular synaptic parameter. The four green symbols (triangle, square, circle, diamond) in the middle panel represent the locations (in terms of parameter values) of the four simulated SPAs shown in the top panel. Bottom, The functional form of the f–I curve of model inhibitory interneurons used for this set of simulations.
Dependence of SPAs on f–I curves for p = 4. A–C, Same as Figure 4 but for p = 4. A, Set of simulations using α = 0.2, λ = 30, nreset = 4. B, Set of simulations using α = 1, λ = 0.98, nreset = 4. C, Set of simulations using α = 3.5, λ = 0.69, nreset = 4. Top, middle, and bottom panel descriptions are the same as given for Figure 4. See Table 3 for units of parameters used above.
A–C, Experimental measurement and quantification of the kinks in f–I curves of CA3 fast-spiking interneurons. A, Top, Infrared-differential interference video microscopy image of a CA3 fast-spiking interneuron. Bottom, Fast-spiking phenotype was determined by the high discharge frequency and high adaptation ratio of the cell. B, Example of a kinky f–I curve from the same cell as in A. C, Example of a less kinky f–I curve. B, C, Top, Actual f–I curve measured. Middle, Slope of the f–I curve 
versus I. Bottom, Change in slope of the f–I curve 
versus I. Local minimum (minima) of this function is (are) the locations of the kink(s) (vertical dotted line) of the original f–I curve. The size of the kink is defined as the ratio of the average slope before the kink to after the kink (see middle).
SPAs not very dependent on mean excitation levels, but more dependent on excitatory conductance fluctuation levels. A–C, Simulations represent progressively increasing mean excitation levels on model inhibitory interneurons. A, ge0 = 0.00383, α = 1, p = 2, λ = 0.81, nreset = 4. B, ge0 = 0.00483, α = 1, p = 2, λ = 0.81, nreset = 4, C, ge0 = 0.00583, α = 1, p = 2; λ = 0.81, nreset = 4. A–C, Top, Summary of occurrence of simulated SPAs. Color codes the same as Figure 4. Bottom, Frequencies of simulated SPAs. Different line colors represent different excitatory conductance fluctuation levels (σe) used for simulation. Green, σe = 0.00102. Deep blue, σe = 0.00152. Pink, σe = 0.00202. Light blue, σe = 0.00252. See Materials and Methods for details on the computation of SPA frequencies and Table 3 for units of parameters used above.
Our simulations were divided into different sets, with each set having fixed cellular properties (see Materials and Methods). Each set of simulations consisted of independent simulations over a range of inhibitory synaptic strengths and a range of excitatory conductance fluctuation levels. To explore how the variation in intrinsic properties of constituent interneurons affected SPAs, we changed the curvature of the f–I curve of individual model interneurons when we moved from one set of simulations to another. Given the equations describing an individual interneuron (equation set 1), to increase the curvature of an individual interneuron's f–I curve, one can increase the initial slope of the curve just above the threshold (increasing α and/or decreasing p), while at the same time increasing the refractory period (decreasing λ and/or increasing nreset). In Figure 3A and B, we show f–I curves of different curvatures as obtained for different α, p, λ, and nreset parameters. Figure 3C shows the firing of a synaptically isolated individual model interneuron for a given set of parameters. As expected, the firing of model interneurons changed when synaptically connected with other model interneurons in the network (compare Figs. 3C, 2, third row of panels).
Robust SPAs occur for particular f–I characteristics of individual interneurons
Figure 4A–C illustrates, respectively, output from three sets of network simulations for a range of gsyn (unitary inhibitory conductance) and σe (excitatory conductance fluctuation) values. In Figure 4, we adjusted the curvature of the constituent interneurons' f–I curve by changing λ, nreset, and α across different sets of simulations (p was fixed at 2 for all the three sets of simulations in this figure). Note that the curvature of the f–I curves of model inhibitory interneurons increases from Figure 4A–C (bottom). These three f–I curves are the same f–I curves as depicted in Figure 3Ai–Aiii, respectively. One can distinguish two approximately linear regimes in each of the f–I curves: the first part when firing first starts and a later region after a kink in the f–I curve. The strength of this kink can be expressed by the slope ratio of the first to the later part of the f–I curve (Table 1). For example, the f–I curve used for model inhibitory interneurons in Figure 4C has a strong kink with close to a vertical slope initially and a later region approaching a horizontal slope, resulting in a slope ratio of 4.2.
The middle panels in Figure 4A–C depict the occurrence of simulated SPAs across synaptic parameter regimes for each set of simulations. One can see that the region of occurrence of SPAs forms a somewhat diagonal band across the gsyn − σe plane, with the “hot” region (yellow to red) occurring toward the lower ends of both the gsyn and the σe axes relative to other “colder” regions. There are some changes of this diagonal band as one increases the curvature of the individual model interneurons' f–I curve (Fig. 4A–C). First, the position of the entire band moves even more toward to the lower-left end of the gsyn–σe plane. Second, the colder regions become more dispersed. In other words, there is a wider parameter regime exhibiting SPAs, but at the expense of the SPAs being weaker in each region for Figure 4C. Thus, it is already apparent from the simulations that is it not only the existence of a kink in the f–I curve but a kink of particular characteristics that is important for the presence of SPAs. This will be made clear in the later section of mathematical analyses.
For each set of simulations in Figure 4A–C, four sample simulations, as given by the green symbol, chosen at different “places” of the gsyn–σe plane (labeled by different shapes of green symbols) are shown with their respective LFP analogs in the top panels. In Figure 4, A and B, one can clearly see the transition from high network activity (green circle) to low network activity (green square) as we increase gsyn, with SPAs occurring at middle gsyn values (green triangle). However, this transition is rather “soft” (meaning that the change from high network activity to low network activity is gradual) for a network in which individual model inhibitory interneurons have a strongly kinked f–I curve (Fig. 4C). Overall, there seems to be a competition between the strength of SPAs and the area of the gsyn–σe parameter space for which simulated SPAs occur. A network with model interneurons having a moderately kinky f–I curve (Fig. 4B) seems to offer the best compromise between strength and pervasiveness of SPAs for inhibitory strength and excitatory fluctuation balances.
Since p determines the nature of bifurcation of a neuron at threshold (see Materials and Methods), it is a potentially important quantity affecting SPAs. In general, f–I curves using p = 4 have smaller curvatures than those using p = 2. Three examples of f–I curves using p = 4 are shown in Figure 3B.
In Figure 5, we depict the occurrence of simulated SPAs as in Figure 4, except here we use p = 4. The transform of the diagonal band on the gsyn–σe plane as one increases the curvature of an individual neuron's f–I curve while keeping p = 4 a constant (Fig. 5A–C) follows a somewhat similar pattern to that of case p = 2. However, in terms of both strength of SPAs and extent of occurrence on synaptic parameter regimes, the p = 4 case is weaker than the p = 2 case (although in Fig. 5A one can still notice some strong SPAs occurring in a small stretch of gsyn–σe plane). In fact, it is somewhat difficult to discern SPAs with the naked eye in some of the colder regions for the p = 4 case. We also performed a set of simulations using p = 1 (data not shown) and found that SPAs were barely perceptible for the explored parameter regimes.
Generation of SPAs requires a kink in the f–I curve of interneurons
A major difference between f–I curves with p = 4 relative to p = 2 is the prominence of the kink for curves with p = 2 but not with p = 4. Table 1 shows the calculated kink size in terms of slope ratio and location for each of our six model interneurons. It is clear from Figure 3A,B and Table 1 that the kinks in f–I curves with p = 2 are more prominent relative to those in curves with p = 4. The particular f–I curves of model interneurons in which network simulations produce the most robust SPAs (Fig. 4A,B) have kink slope ratios close to 2 but much less than 4 (Table 1).
Together, our results suggest that threshold behaviors (as manifest by the values of p) of individual neurons are critical in controlling the occurrence of SPAs. Specifically, given the resolution of our simulations, there must be a moderate kink in the f–I curve (i.e., where the slope of the curve decreases more abruptly), with the kink in approximately the gamma frequency regime for SPAs occurring most robustly. The slope ratio of this optimal kink should be ∼2–3.
The idea that a kinky f–I curve of interneurons is instrumental for SPAs is further reinforced when we later apply our analyses using mean-field equations (equation set 6). We note that the much used inhibitory interneuron model of Wang and Buzsáki (1996) has kink characteristics close to these specifications (kink slope ratio of ∼2 at ∼12 Hz). Thus, inhibitory model networks of Wang–Buzsáki-type interneurons can produce SPAs (similar to those of Fig. 4A). However, if the kink becomes too large (as in Fig. 3Aiii), then SPAs again become much less robust (as shown in Fig. 4C).
Effects of network heterogeneity and non-all-to-all coupling on the occurrence of simulated SPAs
Hitherto we have only used 120-cell all-to-all coupled networks of inhibitory interneurons in which the individual interneurons are homogeneous with the same intrinsic characteristics. Clearly, biological inhibitory networks are neither homogeneous nor coupled in an all-to-all fashion. Would our SPA observations break down in the presence of heterogeneity and non-all-to-all coupling? To address this, we performed several additional sets of simulations in which the all-to-all and homogeneous assumptions were relaxed.
We performed >3000 simulations in which heterogeneity was introduced. For these simulations, we introduced a cellular level of heterogeneity in which there was at least a 3% “baseline” heterogeneity in two of the parameters that controlled the f–I curvature (α and nreset) of individual model interneurons. If we now introduced heterogeneity (10–20%) to the previously uniform gsyn value (i.e., still all-to-all coupling, but interneurons might have had different strengths of connection with one another), we found that SPAs still occurred. Our simulations indicated that SPAs tolerated network heterogeneity well, even with 20% heterogeneity in gsyn values. If we now also changed the all-to-all coupling by randomly deleting some connections between interneurons, we found that the occurrence of SPAs was quite sensitive to the level of connectivity. That is, a sparsely connected network did not support SPAs. Specifically, the simulated network required a minimum of 85% connectivity for robust SPAs to occur.
Given our observations regarding the importance of kinky f–I curves, we also performed simulations in which the 120-cell model network had interneurons with two sets of intrinsic characteristics—one (α, p, λ, nreset) set leading to higher kink slope ratios and another set resulting in lower kink slope ratios. Each set had baseline heterogeneity, as mentioned above. We investigated several “pairings” of intrinsic characteristics in different sets of network simulations. For these split-population simulations, we maintained at least 90% connectivity and had at least 10% gsyn heterogeneity. From these simulations, we found that if the interneuronal subpopulation in which the f–I curve had the higher kink slope ratio constituted more than ∼50–75% of the entire population, it became “dominant” (data not shown). That is, the firing activities of the rest of the population were suppressed, and the dominant subpopulation behaved as if it were a single population. The firing patterns of the dominant subpopulation resembled that of the simulations of a single population in Figures 4 and 5. The degree of dominance depended on the difference in the slope ratios of the kinks of the f–I functions of the two interneuronal subpopulations. The greater the difference in the kink slope ratios between the two subpopulations, the higher the dominance (data not shown). Therefore, we found that the occurrence of SPAs was preserved if the majority of the interneurons in the inhibitory network had the “correct” intrinsic characteristics (meaning kink slope ratios between 2 and 3, with the kink located in the gamma range). For the rest of the population, intrinsic characteristics are such that their f–I functions are less kinky (similar to Fig. 3Bi).
We note that although the extent of our heterogeneous simulations is not as exhaustive as that of our homogeneous ones, we think that it is sufficient enough to allow us to conclude that the occurrence of SPAs requires a high level of connectivity and that some, but not necessarily all, interneurons must have kinky f–I curves.
Experiments suggest that individual CA3 fast-spiking interneurons have the required f–I characteristics as predicted by simulations
Given the significant dependence of SPAs on f–I curve characteristics in our simulations, we decided to undertake a set of experiments to determine whether inhibitory interneurons exhibited appropriate f–I curves. The individual model interneurons in our networks are highly simplified single-compartment models that are meant to represent fast-spiking perisoma-inhibiting interneurons in CA3 hippocampus. Moreover, while f–I curves of inhibitory interneurons have been reported in the literature, such as by Hu et al. (2010), the critical aspect in our predictions is that the f–I function of the majority of the fast-spiking interneurons responsible for SPAs should have a kink at frequencies in the gamma range. Such curvatures of f–I functions are not typically reported in the literature, and, indeed, fine current steps are required to allow such an examination.
We performed whole-cell current-clamp recordings from fast-spiking CA3 interneurons (n = 4) and injected 1 s long depolarizing current pulses. The amplitude was stepwise increased by 5 pA to guarantee high f–I resolution. All recorded cells displayed maximal discharge frequencies above 100 Hz and only minimal adaptation (Fig. 6A) (mean adaptation ratio, 1.05). Two of four tested cells exhibited a moderate kink (Fig. 6B) (mean slope ratio, 2.1) in the gamma frequency range (mean frequency of the occurrence of the kink, 55.3 Hz). The other two cells had only weak kinks (slope ratio, <1.5) (Fig. 6C). Additionally, 12 fast-spiking cells from DG and CA1 regions were also analyzed. More than half of them exhibited appropriate kink characteristics, and the rest had only weak kinks (data not shown).
In summary, our data suggest that f–I characteristics of CA3 fast-spiking interneurons form a heterogeneous population. Two of four cells had the kink characteristics, as predicted by our simulations (i.e., location of the kink between ∼25 and 80 Hz; slope ratio of the kink between ∼2 and 3). The remaining two cells had less kinky (i.e., more linear) f–I characteristics, which is consistent with our heterogeneous simulations. Taken together, experimentally determined f–I characteristics of CA3 fast-spiking interneurons support the predictions of underlying mechanisms for the generation of hippocampal SPAs.
Existence and frequency of simulated SPAs are not sensitive to mean excitatory synaptic conductance levels
In the above simulations, we varied the excitatory fluctuation levels within each set of simulations, but the mean excitatory synaptic conductance (ge0) was kept the same. As can be seen (Fig. 4), the level of excitatory fluctuations must be balanced with the level of synaptic inhibition. But what about mean excitation levels? Here, we show results where we varied the mean excitation to the individual inhibitory interneurons of the network. In the simulations described in this section, we reverted to an all-to-all connectivity with uniform gsyn values between cells, and homogeneous intrinsic characteristics of individual interneurons within a simulation. We used the network in which the individual inhibitory interneurons had a “moderate” f–I curve (Figs. 3Aii, 4B), and we decreased or increased ge0 in steps of 0.001 mS/cm2. Our simulations show (Fig. 7) that a stepwise increase or decrease from our standard ge0 value has only a minor influence on the simulated SPAs. For these sets of simulations with varying ge0, one can see that there is not a significant change either in the shape of the “band” representing the occurrence and the strength of SPAs (Fig. 7, top) or in the frequency statistics of SPAs (Fig. 7, bottom) as ge0 is varied. In other words, an increase or a decrease in the overall mean excitation level ge0 by as much as 22% does not drastically change the results of the network simulations. This observation indicates that SPAs are not heavily dependent on the mean level of excitation. However, we note that SPAs are quite dependent on excitatory fluctuation levels (σe). This can be seen in Figure 4 or Figure 7—SPAs can disappear if σe is increased or decreased by only quarter as much (0.00025 mS/cm2) as the changes made to ge0 (i.e., a one-step change defined above). That is, moving vertically up or down by 0.00025 mS/cm2 from a colored region in the gsyn–σe plane (where SPAs are present) moves one into a black region in which no SPAs are present. In other words, simulated SPAs are much more sensitive to changes in excitatory fluctuation levels (σe) than to mean excitatory conductance levels (ge0).
An increase in excitatory fluctuation causes an increase in frequency but a decrease in amplitude of simulated SPAs
If we now focus on the frequencies of simulated SPAs as shown in the bottom panels of Figure 7, we observe the following. For each set of simulations, in general, we see the SPA frequencies decrease as gsyn increases for a given σe (excitatory fluctuation level). However, if we keep gsyn constant while increasing the excitatory synaptic fluctuation level (σe), we will in general see an increase in frequency of SPAs (Fig. 7, bottom panels), despite a general decrease in strength of SPAs. One can clearly see the decrease in SPA strength (as we increase σe) by moving vertically up from any colored region (and seeing that the color gets colder or turns dark altogether) in any of the top panels of Figure 7. The increase in frequencies with increasing excitatory fluctuations is more noticeable at the lower end of excitatory fluctuation level, partially due to the simplistic way SPA frequencies are computed. We also find that there is an increase in frequency as we increase the curvature of the individual inhibitory neurons' f–I curve (data not shown).
Simulations suggest that modulation of SPAs by adenosine antagonists is due to an increase in excitatory fluctuations and not to mean excitatory levels
Our network simulations are motivated and designed due to the in vitro spontaneous SPAs that robustly occur in hippocampal slices (Wu et al., 2005b). Can we apply our simulation results to understand as yet unexplained observations pertaining to hippocampal SPA experiments? Wu et al. (2009) have recently shown that adenosine antagonists affect hippocampal SPAs by increasing their frequency, lowering their amplitude, and increasing their variability (when compared with hippocampal SPA events of the same slice before perfusion of adenosine antagonists). Adenosine stimulates A1 receptors and, hence, presynaptically inhibits glutamatergic synapses. Thus, adenosine antagonists would lead to increased glutamatergic drive. However, it is unknown (from experiments alone) whether this is due to changes in tonic or phasic excitation. Assuming that mechanism(s) underlying our network simulations producing SPAs are in operation in hippocampal networks, we predict the following (given the increase in frequency and the decrease in amplitude of SPAs seen in the simulations with increasing excitatory fluctuations): hippocampal SPAs are modulated by adenosine antagonists due to an increase in excitatory fluctuation levels rather than to mean excitatory drive changes.
Mathematical analyses using mean-field equations show that network multistability underlies the emergence of SPAs
We have shown that SPAs can be produced in model networks of fast-spiking interneurons if f–I curves of the interneurons exhibit the appropriate kinky behavior. Furthermore, we find that in experiments, fast-spiking interneurons do have kinky f–I curves as predicted. But why does changing the shape and functional form of individual inhibitory interneuron f–I curves have such a profound effect on the occurrence of SPAs? In other words, what is the underlying mechanism that allows SPAs to occur? The goal of this section is to show that the underlying mechanism of SPAs is network multistability. Given the coexistence of multiple network firing states (i.e., network multistability), switching of the network between these states due to excitatory fluctuations would bring about the slow timescales to produce SPAs. We show the existence of multistability by approximating our network simulations with a mean-field model. The advantage of the mean-field approximation is that it allows us to simplify the mathematics and therefore better visualize the underlying mechanisms for SPAs. We assume that the network breaks into one firing and one nonfiring cluster of interneurons. The size (number of interneurons) of the firing cluster is Nb, while that of the nonfiring cluster is Na. This assumption allows us to drastically reduce the dimensionality of our system of equations and to predict how the occurrence of stable clusters of firing interneurons depends on different parameters. Further details are provided in Materials and Methods, but the essence boils down to the analyses giving stable solutions of firing states of different Nb values on the gsyn line. Stable solutions are defined as solutions having negative eigenvalues (units of eigenvalues in ms−1).
Our simulation results in Figures 4 and 5 show colored plots on (gsyn – σe) planes, indicating that SPAs emerge when synaptic inhibitory conductances, gsyn, and excitatory fluctuation balances are appropriate, given that f–I interneuron curves are moderately kinky. We performed mean-field analyses on these model simulations, and Figure 8 illustrates the results of our analyses. The analyses are shown in Figure 8 on (gsyn − log2 Nb) planes, to allow a correspondence with gsyn, which are the x-axes of Figures 4 and 5. We note that because the mean-field model approximates the simulations, the gsyn values for which SPAs occur in the simulations are not expected to be identical to the gsyn values for which multistability occurs in the analyses. Figure 8A–C depicts the results of our mean-field analyses for the simulations in Figure 5A–C, respectively (f–I curves for p = 4), while Figure 8D–F depicts mean-field analysis results for simulations in Figure 4A–C, respectively (f–I curves for p = 2). For convenience, the corresponding f–I curve is included as an inset on the lower right side of each plot in Figure 8A–F.
Mean-field model analyses explain the generation of SPAs. A–F, Mean-field solutions that approximate the full simulations are shown on (gsyn − log2 Nb) planes, where Nb is the number of firing cells in a cluster. The lines represent where stable solutions exist (i.e., eigenvalues being negative) for given gsyn values shown on the x-axis. For Nb = 1, the largest eigenvalue is negative as soon as there is a solution, and it remains constant throughout. Therefore, the line for Nb = 1 extends all the way to infinity. The inset shows the f–I curve used in network simulations in which the mean-field equations are approximating. The pair of dotted lines (D) illustrates an example of the overlapping of two different stable network firing states. Mean-field solutions corresponding to simulations with f–I curve given by the following: (A) Figure 3Bi: α = 0.2, p = 4, λ = 30, nreset = 4; (B) Figure 3Bii: α = 1, p = 4, λ = 0.98, nreset = 4; (C) Figure 3Biii: α = 3.5, p = 4, λ = 0.69, nreset = 4; (D) Figure 3Ai: α = 0.41, p = 2, λ = 20, nreset = 4; (E) Figure 3Aii: α = 1, p = 2, λ = 0.81, nreset = 4; (F) Figure 3Aiii: α = 10, p = 2, λ = 0.74, nreset = 32. See Table 3 for units of parameters used above.
To understand the relationship between SPAs and individual cellular characteristics, we investigated how the relative distribution of mean-field solutions of stable (i.e., negative eigenvalue) firing states was affected by the f–I curves of individual interneurons. By distribution, we mean the distribution of gsyn intervals in which different solutions to the mean-field model are stable. The lines on each of the six plots in Figure 8 denote the gsyn intervals on which the mean-field solution (equation set 6) for a particular number, Nb, of firing cells (y-axis is log2 Nb) is stable. If more than one line is present for the same gsyn values (i.e., there is an overlap), then there exists more than one possible stable cluster of firing interneurons for the same gsyn value. In other words, multistability exists. We note that successive log2 Nb values have gsyn intervals that “stack” in a “diagonal” manner, with smaller Nb values occupying larger gsyn values. Larger gsyn values generally lead to mean-field solutions with a smaller number of firing interneurons (i.e., smaller Nb). Intuitively, this is because a larger gsyn gives rise to stronger inhibition so that interneurons could more easily suppress other interneurons in the networks. Moreover, if there are fewer interneurons firing, each firing interneuron has to fire more frequently to suppress other interneurons that are not firing. This was observed in our simulations in which, during the low activity state, each firing interneuron generally fired at a higher frequency than interneurons firing during the high activity state. Therefore, one can generally associate mean-field solutions of smaller Nb values with higher firing frequencies of individual interneurons.
We now go through the mean-field results of Figure 8A–F, starting with Figure 8A. When f–I curves of individual model interneurons are almost linear (Fig. 8A), mean-field solutions for different Nb values have similar gsyn interval lengths (except for Nb = 1). This can be intuitively understood as follows: since the derivative of the f–I curve is associated with the stability of mean-field solutions (see Materials and Methods), an overall more linear f–I curve confers similar degrees of stability for mean-field solutions of different Nb values (because the derivative, df/dI, is approximately constant over the f–I curve). Thus, similar gsyn intervals or lengths of lines can be seen in Figure 8A. Moreover, because there is little curvature (i.e., second derivative, d2f/dI2, is close to zero) of the f–I curve in Figure 8A, the gsyn intervals of different Nb values overlap moderately so that network multistability exists and we expect to observe SPAs. In this case, the overlap exists for both large and small Nb values. As such, SPAs are expected to be “strong” (i.e., of larger magnitude) because it was possible for the network to switch to a state in which many interneurons (high Nb values) are firing. However, SPAs are expected to occur over a relatively small stretch of gsyn values because the overlap is concentrated over a small range of gsyn values. This is the case (i.e., stronger SPAs over a narrower range of gsyn values) as we can see in the network simulations of Figure 5A in which individual interneurons had f–I curves of Figure 8A or Figure 3Bi.
Going from Figure 8A to 8B,C, we increase the curvature of f–I curves of individual interneurons. We observe that mean-field solutions with fewer firing interneurons (smaller Nb) occupy a larger interval of gsyn values (Fig. 8B,C), while those with more firing interneurons (larger Nb) occupy smaller intervals. There is also less overlapping in gsyn intervals of the mean-field solutions, and the diagonal “stacking” of the mean-field solutions on the gsyn intervals is more slanted and bent away from the vertical (Fig. 8B,C). To understand this, note that there is a larger curvature (i.e., non-zero |d2f/dI2| values) along the f–I curve in Figure 8, B and C; thus, mean-field solutions for different Nb values no longer have similar gsyn intervals in Figure 8B,C as they do for the f–I curve in Figure 8A. In particular, we note that with the initial portion of the f–I curves (Fig. 8B,C) being moderately steeper than the later portion, small to intermediate Nb values have increased gsyn intervals, whereas the intervals for larger Nb values are decreased. Furthermore, the relative constant curvature of f–I curves with p = 4 as compared to p = 2 (i.e., smaller variation of d2f/dI2 over the f–I curve for p = 4 relative to p = 2) ensures little overlap of mean-field solutions in gsyn values (mean-field solutions of different Nb values are uniformly separated from one another on the gsyn line). Thus, it is expected that SPAs would be difficult to observe using f–I curves shown in Figure 8, B and C. This is the case as seen in Figure 5, B and C.
Figure 8D–F shows mean-field results for f–I curves with p = 2 (p is the variable that determines the characteristics of the interneuron f–I curve just above threshold; see equation set 2). There are two differences between f–I curves with p = 4 (Fig. 8A–C) and those with p = 2 (Fig. 8D–F). First, with p = 2, the f–I curves are generally of a higher curvature (i.e., higher |d2f/dI2| values) than with p = 4. Second, the variation of curvature of f–I curves with p = 2 occurs over a narrower range of injected current values, as opposed to f–I curves with p = 4, where the curvature varies more evenly over a wider range of injected current values. As a result of these two characteristics, a kink is more prominent for f–I curves with p = 2 but less so with p = 4. One of the effects of the kink in the f–I curve is that it introduces an uneven separation in gsyn intervals for different Nb values. For small to intermediate Nb values (Nb ≤ 8) for which interneurons fire at a higher frequency, the gsyn intervals increase in length, allowing more overlap (Fig. 8D,E) because the portion of f–I curve above the kink is substantially less steep than the portion below the kink. This is where one finds robust coexistence of stable firing states of low to intermediate Nb values responsible for SPAs. Further increasing the kink of individual interneuron f–I curves (Fig. 8F) nudges solutions with intermediate Nb values into the periphery; thus, only solutions with small Nb values dominate the gsyn line. As a result, we see SPAs of low strength across a wider range of gsyn values. This is shown in Figure 4C, which uses the f–I curve depicted in Figure 8F or Figure 3Aiii.
To further illustrate the relationship between SPAs and network multistability, we performed a quantification of the amount of overlap. In the last column of Table 1, we show the amount of overlap in gsyn values for mean-field solutions with small to intermediate Nb values for each of our six mean-field analyses (Fig. 8A–F). The overlap is defined as the total length in units of gsyn in which there is a coexistence of any two stable firing states of Nb ≤ 8. It is clear that the gsyn overlap is maximized for f–I curves of p = 2 and with the kink slope ratio between ∼2 and 3. These f–I curve characteristics, when used to represent individual interneurons in network simulations, produced the most robust SPAs (Fig. 4A,B). There is also considerable overlap of gsyn intervals in which an almost linear f–I curve was used (Fig. 8A). In the corresponding set of network simulations (Fig. 5A), we see that SPAs of high strength occurred over a narrower range of gsyn values. For the three mean-field solutions with the smallest overlap (Fig. 8B,C,F), we see only minimal SPAs in the corresponding sets of simulations (Figs. 5B,C, 4C, respectively).
In summary, from the above analyses, we found that the occurrence of SPAs in network simulations was closely associated with the degree of coexistence (overlap) in mean-field solutions with low to intermediate Nb values. As shown in our simulations, movement between different states occurred in the context of appropriate levels of excitatory fluctuations and inhibitory strengths. Specifically, the two f–I curves in which there was substantial overlap of low to intermediate firing states for the mean-field solutions coincided with the two sets of full simulations in which most robust SPAs were observed (Table 1). Moreover, the f–I characteristics (i.e., kink size and location) of model interneurons in which robust SPAs were demonstrated in simulations were also shown to exist experimentally in interneurons of CA3 hippocampus (Fig. 6). Although the mean-field analyses only approximate the network simulations, the correspondence between the analyses and the simulations as described above is robust. As such, the requirement of network multistability as an underlying mechanism for the existence of hippocampal SPAs is shown.
To make our explanations more concrete and to show further correspondence with the experiment, we bring forth several supporting points and present an idealized figure. Figure 9A depicts an idealized arrangement of stable states to illustrate inhibitory network multistability, which occurs when f–I curves are moderately kinky. As shown, there are stable states with smaller and larger numbers of firing cells, Nb, for the same gsyn values. Thus, at a given time, the network may reside in any stable state. In our proposed mechanism, SPAs are due to movement between stable intermediate network activities and stable low network activities. In experiments on hippocampal SPAs, we observed low and intermediate firing states. That is, some inhibitory interneurons fired between SPA episodes in CA3 hippocampus, and many more interneurons fired during SPA episodes (Wu et al., 2005b, their Fig. 7B). Furthermore, this can also be clearly seen in our model simulations; see raster plots of Figure 2 in which model interneurons fired between SPA episodes.
A, B, Schematics showing idealized arrangement of mean-field solutions for the occurrence of robust and pervasive SPAs, and f–I curve showing required locations for gamma output and SPAs. A, Schematic illustrating an idealized arrangement of stable cluster solutions for the robust occurrence of strong SPAs—multiple coexistence of stable states that include low to intermediate Nb values—as illustrated by the region between the two dotted lines. B, f–I Curve illustrating difference between gamma and SPAs. For SPAs, the network essentially operates at the kink region. For gammas, the network operates at the latter more linear region of the f–I curve, making it less likely that movement between low and high firing states can occur.
Our mechanism shows the requirement for suitably kinked f–I curves at gamma frequencies. This means that during intermediate network activities when many cells were firing, they should have been doing so at gamma frequencies since they had to be located near the f–I kink to allow movement between states due to excitatory fluctuations. This required “location” is illustrated by a dashed rectangle on an appropriate kinky f–I curve in Figure 9B. We examined several simulated SPA events. We looked at ∼50 SPA events from two different simulations (that had SPAs similar to those of our experiment) and found that for each SPA event (each of ∼200 ms duration), approximately half the model interneurons fired only one spike and the rest fired two or more spikes. Of the multiple firing spikes during SPAs, the average frequency was ∼45 Hz (i.e., in the gamma range). In our experiment, a mean frequency of ∼70 Hz is reported (see Wu et al., 2005b, their Fig. 7). In other words, the spiking model interneurons during SPAs sit near the kink of the f–I curve so that they are able to move between firing or not firing with appropriate excitatory fluctuation (noise) levels.
Comparison between gamma and hippocampal SPAs
We have described how SPAs can arise in inhibitory networks via network multistability. However, inhibitory networks are also known to exhibit gamma oscillations that are of much higher frequencies (25–100 Hz). How is it that the same inhibitory network can sustain oscillations of apparently two different ranges of frequencies? The answer lies in the different inhibitory conductance values and excitatory fluctuation levels required by each of gamma and hippocampal SPAs. Roughly speaking, one can think of SPAs as “another kind” of gamma, because whenever a collection of interneurons in the network fires during SPAs, each neuron fires at approximately the gamma frequency. The difference between SPAs and gamma is that in gamma, the number of interneurons participating in firing is stable. In SPAs, this number is modulated by an emergent timescale. That is, at times there are more interneurons firing and at other times there are fewer interneurons firing, thus creating a dominant lower frequency (frequency of modulation of gamma) in the power spectrum. For SPAs to occur in the first place, inhibitory relative to excitatory input has to be such that the network, on balance, essentially operates at the kink part of the f–I curve (Fig. 9B). The excitatory fluctuations then bring about the SPAs by allowing movement between different stable network states. On the contrary, for gamma oscillations generated by the interneuron network gamma mechanism (Whittington et al., 2000; Tiesinga and Sejnowski, 2009), one in general requires smaller inhibitory conductances and lower excitatory fluctuation levels than are required by SPAs. The interneuronal network during gamma is more likely to behave as a single cluster, and there is a stronger synchronization of spike timings between individual interneurons during gamma oscillations. Individual interneurons essentially fire at frequencies located at the farther and less steep end of the f–I curve beyond threshold (Fig. 9B, dashed circle). Since the excitatory fluctuations are not strong enough to move the system out of a firing state, it is difficult for slow timescales to emerge.
Discussion
Summary of results
We have used a combination of experiments, data extraction, and mathematical modeling techniques in previous (Wu et al., 2005b; Ho et al., 2009) and present work to elucidate the mechanisms underlying a form of SPAs. These SPAs originate in CA3 hippocampus and are due to coherent inhibitory activities. Our findings show that hippocampal SPAs are an emergent phenomenon, and the intrinsic properties of individual interneurons are salient determinants of the phenomenon. We summarize our hippocampal SPA mechanism as follows: a suitably kinky f–I curve of interneurons promotes maximal coexistence of network states between low and intermediate numbers of firing cells (Fig. 9A, region between the pair of dotted lines). SPAs occur when the excitatory fluctuations are strong enough to move the network from one firing state to another, thus creating the population slow timescales we observed in hippocampal SPA experiments. The properties at the initial portion of individual interneuron f–I curves control most characteristics of hippocampal SPAs. For hippocampal SPAs to exist, the f–I curves of interneurons require a kink at the gamma frequency range above threshold (Fig. 9B, dashed rectangle). We have also experimentally measured the f–I characteristics of CA3 fast-spiking interneurons and quantified their kinks, finding that they are consistent with our model predictions.
While the concept of network multistability in neuroscience is not novel (van Ooyen et al., 1992; Freyer et al., 2009; Fröhlich et al., 2010), to the best of our knowledge we are the first to have elucidated the conditions that promote network multistability. From our simulations we have demonstrated that SPAs are more sensitive to excitatory conductance fluctuations than to changes in mean excitatory levels. We have applied this result to understand the effects of adenosine antagonists on hippocampal SPAs.
Our proposed model of hippocampal SPAs is based on simulations and analyses independent of specific ionic mechanisms responsible for the kinkiness of f–I curves of interneurons. Therefore, any ionic mechanism that produces the correct curvature of f–I curves can support SPAs in theory. However, since the kink of the f–I curve is primarily a result of the rapid onset of spiking of the interneuron just above threshold and the refractoriness of the interneuron, ionic mechanisms responsible for these two characteristics are probable candidates for sustaining SPAs. The fast-spiking aspect produces a steep segment at the very initial portion of the f–I curve just above threshold. This leads to a reduction in the formation of large firing clusters (large Nb). The refractoriness ensures a less steep segment at the ensuing portion of the f–I curve, which promotes overlapping of firing states with low to intermediate Nb values. For example, voltage-gated potassium channels Kv3.1 and Kv3.2, which are implicated in controlling the fast-firing characteristics of inhibitory interneurons (Lau et al., 2000; Rudy and McBain, 2001; Vyazovskiy et al., 2002) would be candidates for modulating the occurrence of SPAs.
Hippocampal SPAs as meta-organization of inhibitory networks
A recurring theme in many natural phenomena is how simple elements coupled with simple rules of interaction can give rise to complex large-scale structures (Wolfram, 2002). In many cases, systems achieve complexity through self-similarity (Bak, 1996) and its associated “power-law” behavior. There are many examples in nature and in the brain (Mandelbrot, 1984; Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003; Schneidman et al., 2006; Tang et al., 2008; Bullmore et al., 2009; Shew et al., 2009; He et al., 2010; Werner, 2010) in which power-law behaviors have been discovered. These observations signify the organization of individual neurons to represent themselves at a scale beyond what is typically expected of any single neuron. However, physiological mechanisms through which individual neurons self-organize are poorly understood. Currently, a common answer is balanced excitation and inhibition (Shew et al., 2009). However, no clear context is given about this balance. Our work provides a context for the case of inhibitory networks. Our results suggest that, as well as synaptic properties, intrinsic properties of inhibitory interneurons can play a vital role in bringing the network into a state of phase transition. In our case, the self-organization is built upon a “first-order”-like phase transition in which several stable states overlap on a range of gsyn values. The slow timescales are the result of the system “jumping” back-and-forth between these states (Fig. 10A). The self-similarity is imperfect because of the limited number of overlapping states. In a genuine self-similar system, the phase transition is “second-order” (Landau and Lifshitz, 1980; Pathria, 2001; Plischke and Bergersen, 2006) with a single point (e.g., a single gsyn value) onto which every state converges. However, this is not supported by experimental data, because one would then need all the interneuron f–I functions to be more linear and similar to Figures 3Bi or 6C for the network to experience a second-order-like phase transition. It is also biologically infeasible unless extra mechanisms (Bak et al., 1987; Bak, 1996; Levina et al., 2007, 2009; Millman et al., 2010) are in place to overcome the fine-tuning problem. Figure 10B shows the hypothetical arrangement of mean-field solutions corresponding to a second-order-like phase transition when every interneuron in the network has an almost perfectly linear f–I curve above threshold.
A–C, Power law behavior of model SPAs and the coexistence of hippocampal SPAs and SPWs in some hippocampal slices (A). Power spectrum of one of the SPA simulations on log-log scale. Parameters used from simulation: p = 2, α = 0.41, λ = 20, gsyn = 0.041, σe = 0.00102, nreset = 4. Note the approximate linearity of the spectrum for >3 decades (from 0.01 to ∼35 Hz). See Table 3 for units of parameters used. B, Hypothetical arrangement of mean-field solutions of a second-order-like phase transition. This arrangement of mean-field solutions happens when every constituent interneuron in the network has an almost perfectly linear f–I curve above threshold. C, LFP recording on experiments in which there was coexistence of hippocampal SPAs and SPWs. During the LFP recording, the SPW episodes were usually of much larger amplitude (up to a few mV) and much less frequent than hippocampal SPA episodes. Three SPW episodes in this recording are each labeled by a diamond, and two examples (out of many) of hippocampal SPA episodes are marked by crosses. The SPW events in this recording are artificially truncated at 1 mV in this figure.
Extension of our results and related works
We discuss two examples beyond our current experimental setup for hippocampal SPAs in which our findings can shed light. The first one concerns experiments by Wu et al. (2005a). In some hippocampal slices prepared by this group, hippocampal SPAs were found to coexist with another kind of population activity called sharp waves (SPWs). SPWs are the result of collective firing of excitatory pyramidal neurons. Figure 10C depicts a LFP recording in which SPWs coexist with hippocampal SPAs. While there is no clear trend in terms of frequency and amplitude of successive hippocampal SPA episodes in hippocampal slices without the existence of SPWs, here the hippocampal SPA episodes increase in frequency but decrease in amplitude as the next SPW episode is imminent (Fig. 10C). If our results about the effects of excitatory fluctuations on hippocampal SPAs can be carried over to this situation, the behaviors of hippocampal SPAs during two SPW episodes might reveal temporally increasing excitatory fluctuation levels as the hippocampal network approaches the next SPW episode. Since increased excitatory fluctuation levels are indicative of the imminent onset of SPW episodes (de la Prida et al., 2006), our results provide a parsimonious explanation that naturally links the two observations (SPWs and hippocampal SPAs) together (Ellender et al., 2010 described the possible roles of GABA(A) receptor-mediated inhibition in SPWs). Similar ideas might also be used to explain the change in properties in SPAs in hippocampal slices of a Rett syndrome mouse model (Zhang et al., 2008).
During neocortical slow oscillations (Steriade et al., 1993; Timofeev et al., 2000; see Wolansky et al., 2006 for slow hippocampal oscillations), the neocortex sustains a network UP and a network DOWN state. The DOWN state is characterized by the paucity of synaptic activities in the network, while the UP state is the opposite. The network switches between UP and DOWN state with a frequency <1 Hz. It is not fully understood how the neocortex can transition between the two states. There are, however, computational models that reproduce some of the behaviors observed in experiments (Compte et al., 2003; Holcman and Tsodyks, 2006). Although we have not explicitly included pyramidal cells into our simulation, our work raises the possibility that network multistability can be built up from individual neurons with only simple characteristics. One can envision that a minimal computational model of UP and DOWN states can be constructed by adjusting the f–I curves of some mathematical representation of neurons. Network stable UP states would have firing frequencies of individual neurons located in the “not-so-steep” regions of f–I curves. Noises due to synapses or heterogeneity of neuronal population might contribute to the transition between the UP and DOWN states.
Limitations and future work
Although we have included the essential physiological ingredients for the emergence of SPAs, our model has limitations. Besides the simplified single-compartment model representation of individual interneurons, most notable is the absence of electrical synapse representations (gap junctions) between interneurons. Experiments and computational studies have shown that electrical and chemical synapses work together to affect network output (Skinner et al., 1999; Tamás et al., 2000; Traub et al., 2001). Since electrical synapses are believed to enhance the synchronization of spike timing, we suspect that larger excitatory fluctuation levels might be required to sustain SPAs were gap junctional effects included in our simulations. We do not believe our conclusions would be significantly affected. We plan to incorporate excitatory cells directly, expand the inhibitory network size to include more inhibitory cell types, and examine more realistic types of network connectivity architecture.
We note that even though excitatory populations are only implicitly included in our present models so that feedback inhibition is absent, we do not expect drastic changes regarding SPAs if excitatory cells are explicitly included. This is because SPAs are relatively insensitive to mean excitatory drive (Fig. 7) and this is what would be primarily modulated with the inclusion of inhibitory feedback.
We have developed a mean-field theory for our mathematical analysis of the emergence of SPAs. One deficiency of mean-field theory is the incorrect representation of fluctuations. The development of an “effective potential” (similar to free energy in thermodynamics) via a field-theoretical description is one such approach to address the problem (Zinn-Justin, 2002). Buice and Cowan (2007) and Buice et al. (2010) have used this approach in a neural network context. While we believe our results would not be significantly affected, it would be interesting to apply the above approach to our interneuronal networks.
Footnotes
This work was supported by a Natural Sciences and Engineering Research Council of Canada Discovery Grant (L.Z., F.K.S.) and a PGS-D award (E.C.Y.H.); Queen Elizabeth II Graduate Scholarship in Science and Technology of Ontario (E.C.Y.H.); University of Toronto Open Fellowships (E.C.Y.H.); Lichtenberg Award, VW-Foundation (M.B.); BIOSS Centre for Biological Signalling Studies (M.B.); and Excellence Initiative of the German Federal and State Governments (GSC-4, Spemann Graduate School) (M.S.). Computations were performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund Research Excellence; and the University of Toronto. We thank Dr. Nathan Insel for critical reading of the manuscript and the reviewers for constructive comments.
- Correspondence should be addressed to Frances K. Skinner, Toronto Western Research Institute, University Health Network, MP13-317, 399 Bathurst Street, Toronto, ON M5T 2S8, Canada. frances.skinner{at}utoronto.ca
