Abstract
Computational studies have suggested that stochastic, deterministic, and mixed processes all could be possible determinants of spontaneous, synchronous network bursts. In the present study, using multicellular calcium imaging coupled with fast confocal microscopy, we describe neuronal behavior underlying spontaneous network bursts in developing rat and mouse hippocampal area CA3 networks. Two primary burst types were studied: giant depolarizing potentials (GDPs) and spontaneous interictal bursts recorded in bicuculline, a GABAA receptor antagonist. Analysis of the simultaneous behavior of multiple CA3 neurons during synchronous GDPs revealed a repeatable activation order from burst to burst. This was validated using several statistical methods, including high Kendall's coefficient of concordance values for firing order during GDPs, high Pearson's correlations of cellular activation times between burst pairs, and latent class analysis, which revealed a population of 5–6% of CA3 neurons reliably firing very early during GDPs. In contrast, neuronal firing order during interictal bursts appeared homogeneous, with no particular cells repeatedly leading or lagging during these synchronous events. We conclude that GDPs activate via a deterministic mechanism, with distinct, repeatable roles for subsets of neurons during burst generation, while interictal bursts appear to be stochastic events with cells assuming interchangeable roles in the generation of these events.
Introduction
Characterization of the simultaneous spatiotemporal firing patterns of multiple neurons is fundamentally important in the understanding of both normal and pathological processes within circuits. Synchronous or near-synchronous firing of populations of neurons contributes to coding of activation strength as an input mechanism. Similarly, precise, ordered temporal firing patterns of multiple neurons is important as a delayed input activation mechanism (Szatmáry and Izhikevich, 2010), encoding input strength and possibly contributing to induction of long-term potentiation. The temporal firing pattern itself could be an additional method of information representation within a neuronal circuit, contributing to rapid temporal coding (VanRullen et al., 2005). Using multicellular calcium imaging, activation patterns of multiple cells (>100 neurons) can be studied and analyzed (Mammano et al., 1999; Garaschuk et al., 2000; Stosiek et al., 2003; Hirase et al., 2004; Ikegaya et al., 2004; Cossart et al., 2005). This provides information on functional connectivity and insights into circuit dynamics (Takahashi et al., 2007; Bonifazi et al., 2009; Sasaki et al., 2009; Rothschild et al., 2010; Hofer et al., 2011). Most calcium imaging studies reported in the literature have focused on relatively random, dispersed firing networks as this technique provides a unique source of information on functional connectivity within neuronal microcircuits.
Less attention has been paid to spontaneous or evoked synchronous bursts within networks. In these studies, because of the near total activation of the neuronal circuit, all cells appear “connected,” and little insight can be provided in the elucidation of functional connectivity patterns and input/output mechanisms. However, scrutinizing the cellular activation patterns within synchronous network bursts is critical in understanding burst generation mechanisms and network dynamics (Fellin et al., 2004; Crépel et al., 2007; Allène et al., 2008). Spontaneous synchronous network bursts occur in immature brain during development, and in adult brain as a component of sensory responses. They also are prevalent in pathological conditions such as in epilepsy. In an epileptic context, there has been a significant focus on the concept of pathological plasticity inducing the development of “burster” neurons, which drive epileptic activity in focal regions in a deterministic manner (Jensen and Yaari, 1997; Sanabria et al., 2001; Becker et al., 2008; Chen et al., 2011). If epileptic network bursts are initiated by bursters, identifying those cells is critical in understanding mechanisms of pathologic burst generation. Network bursts could also be generated in a homogeneous population through a stochastic process (Miles and Wong, 1983). Network modeling studies have demonstrated that randomly connected populations of homogeneous neurons can generate synchronous network bursts (Izhikevich, 2003; Bogaard et al., 2009).
In the present study, we investigated cellular firing patterns underlying spontaneous synchronous network bursts in the developing hippocampus by using multicellular calcium imaging and fast confocal microscopy (up to 350 Hz recording speed).
Two types of synchronous network activity, giant depolarizing potentials (GDPs) and interictal bursts, were investigated. We found that GDPs exhibited a nonrandom, deterministic burst structure, with ∼5% of cells being characterized as an “early responder” population with high confidence. In contrast, we found that spontaneous interictal bursts occurring in the presence of bicuculline did not exhibit strong deterministic structure, with no early responder population evident. Bicuculline-induced interictal bursts appeared to be generated in a homogeneous population through a stochastic process.
Materials and Methods
Slice preparation and staining.
Brain slices were prepared as described previously (Carlson and Coulter, 2008). Briefly, the brain of male and female rats or mice [postnatal day 4–9 (P4–P9)] were removed, glued to an agar block, and sectioned (350–400 μm thickness) in ice-cold artificial CSF (ACSF) containing the following (in mm): 130 sucrose, 3 KCl, 1.25 NaHPO4, 1 MgCl2, 2 CaCl2, 26 NaHCO3, and 10 glucose. A dye solution of 10 μl of 0.1% Oregon Green BAPTA-1 AM (Invitrogen) in DMSO, 2 μl of 5% Cremophor EL (Invitrogen) in DMSO, and 2 μl of 10% pluronic acid (Invitrogen) in DMSO in 2 ml of ACSF [containing the following (in mm): 130 NaCl, 3 KCl, 1.25 NaHPO4, 1 MgCl2, 2 CaCl2, 26 NaHCO3, and 10 glucose] was prepared as described in the literature (Takahashi et al., 2007). The slices were submerged in a dish containing 2 ml of a dye solution for 45 min at 37°C in a chamber under humidified 95% oxygen/5% carbon dioxide. Slices were washed and incubated in oxygenated ACSF at room temperature for at least 30 min before imaging.
Imaging.
Slices were transferred to an imaging chamber and perfused with modified ACSF containing an elevated level of potassium (6 mm) during imaging. A Live Scan Swept Field Confocal Microscope (Nikon Instruments) equipped with an Ar/Kr ion laser (Innova 70C; Coherent) was operated from NIS-Elements software (Nikon). Images were captured with water-immersion 40× [numerical aperture (NA) = 0.8] or 16× lenses (NA = 0.80) and one of two EM-CCD cameras (Cascade 512B or Cascade128+; Photometrics). We focused on a region of area CA3 (area CA3b) that generates more spontaneous synchronized events than other areas within the slice (e.g., the dentate gyrus or CA1). The Cascade 512B camera provided acquisition rates of up to 25 Hz with a full frame (512 × 512 pixels), and 89 Hz with 4 × 4 binning (128 × 128 pixels). The Cascade128+ camera provided an acquisition rate of 350 Hz with full frame (128 × 128 pixels) resolution when combined with Swept Field microscopy.
Image analysis.
All time-series movies (ND2 file; Nikon Instruments) were converted into TIF format by an ImageJ plugin (ND2 reader; Nikon Instruments), and regions of interest (ROIs) corresponding to cell somata were selected manually from a maximum response image using ImageJ software. The average intensity of ROIs was calculated and exported as a text file. Since Oregon Green calcium transients are reflected by an increase in fluorescence intensity, we chose to determine ROIs on a maximum response image so as to ensure that the majority of active neurons were included in our analysis. Custom-written Matlab (MathWorks) codes were used to calculate background subtracted normalized fluorescence change (ΔF/F0). Raster plots were then generated using ImageJ software. Onset time detection was performed in Matlab as described below. Color-coded onset time mapping was generated by assigning colored circles in the center of ROIs by Matlab and ImageJ.
Onset time estimation.
The rising phase of synchronous ΔF/F0 traces was sigmoidal in shape. To estimate onset time, this was fitted to a four-parameter sigmoidal equation (Zwietering et al., 1990), where A is the plateau, B is maximum slope, C is a lag time, D is an offset, and e is Euler's constant:
We estimated “onset time” as the time when f(t) reached to 5% of the plateau. To evaluate the estimation error in this curve-fitting approach, we conducted a simulation by adding normally distributed noise with various noise levels. When the noise level (SD of baseline before the response) was 20% of the response signal amplitude, the SD of onset time estimation was 13 ms at a sampling rate of 100 Hz (n = 500). The SD of onset time estimation reduced to 8.8 and 5.4 ms as the noise level was reduced to 15% and 10%, respectively. Note that the average noise levels of our data from GDP experiments were ∼10%, ranging from 5% to 20%. In a subset of experiments in which we combined whole-cell patch-clamp recordings with imaging of calcium transients at our usual temporal resolution, we found that calcium transients occurred within 1 ms of action potential (AP) firing on average, with an SD of 7 ms. For the data obtained at 350 Hz, we used the time frame when ΔF/F0 exceeded the 5% threshold.
Statistical analysis.
Statistical analyses were performed using R2.13 (R foundation for Statistical Computing, http://www.R-project.org) and GraphPad Prism 4. Initial observations suggested that activation time was a function of image pixel location (i.e., bursts exhibited directional propagation). To characterize this propagation, the data were processed by fitting a regression spline with up to 5 df using library (spline) in R. To characterize the degree of correlation among onset time for repeated bursts (waves), we estimated Pearson's correlation coefficient in a pairwise fashion using onset time and/or onset time residuals from the spline fit. We also examined the nonparametric ranking order statistic, Kendall's coefficient of concordance (W), to compare consistency of activation order among multiple bursts for cells that are in close spatial proximity. Kendall's W value is an indicator of rank consistency in cellular activation time among the bursts, and varies from 0 to 1, where W = 1 signifies totally consistent activation order among repeated bursts, and W = 0 is completely inconsistent activation ranks.
To explore the clustering of neuronal response times, or to determine and characterize a group of neurons that represent a cluster of early responders or “followers,” we conducted latent class (LC) model analysis. LC analysis, is a parametric statistical method for clustering data into a number of unobserved, or latent, classes (McLachlan et al., 2002). We know of no formal hypothesis testing procedure in the statistical sciences that can be used to test for the existence of multiple distributions or classes against the null hypothesis of a single underlying distribution (McLachlan and Peel, 2000). However, in the spirit of hypothesis testing, we wished to ensure with high probability that the LC analysis was unlikely to declare more than one class when only one was present. Standard latent class modeling involves choosing the number of classes based on the minimum of information criteria such as the Bayesian Information Criteria (BIC), or in our case an adjusted version of BIC (adjBIC). To characterize the LC analysis procedure, we conducted a simulation experiment based on data, generated with a structure similar to our data, but from either one or two distributions. Then we determined how frequently we chose models with one, two, or three classes. We were specifically interested in documenting how often we chose more than one class when the data had one class. This analysis confirmed that the procedure infrequently (typically no more often than 5% of the time) yielded spurious classes when none existed. The LC model estimates a mean and SD for the observations for each class and assigns a posterior probability of membership in each class to each cell. Cells were assigned to class based on their maximum posterior probability, and we also performed a sensitivity analysis for these results based on a maximum posterior probability in excess of 0.90. LC analysis differs from many other clustering methods in reporting uncertainty of class assignment through the posterior probability. The LC analysis was implemented using library (lcmm) in R.
Results
Calcium imaging and neuronal firing
We performed multicellular calcium imaging of hippocampal CA3 pyramidal neurons with a 40× objective lens, followed by off-line image analysis. With the 40× lens, the image size was 204 × 204 μm, typically containing 50–100 pyramidal neurons within the focal plane of an image. The increase in fluorescence signal in neuronal cell bodies in multicellular calcium imaging is a result of calcium influx accompanying AP firing (Stosiek et al., 2003; Ikegaya et al., 2004; Cossart et al., 2005; Sasaki et al., 2008; Vogelstein et al., 2009). To ensure the correlation between fluorescence signal change and neuronal firing, a patch pipette electrode was positioned in close proximity to a cell and loose-patch recording was performed during calcium imaging. Figure 1A shows a trace of fluorescence change measured from a cell that exhibited spontaneous activity. Figure 1B shows a simultaneously recorded trace of loose patch recording measured in the same cell as Figure 1A. By comparing Figure 1A and Figure 1B, it is evident that the fluorescence and current changes were concurrent, suggesting that fluorescence changes were directly related to AP firing. In additional combined whole-cell patch-clamp recording and cellular imaging experiments at 100 Hz temporal resolution, and estimating calcium signal onset using the curve-fitting approach described in Materials and Methods, we found that the onset of the calcium signal relative to the occurrence of the AP was nearly simultaneous under our imaging conditions (mean difference = 1 ms, SD = 7.4 ms) (Fig. 1C). In addition, loose patch recording reported the number of APs elicited during a burst-firing event. Figure 1D plots the correlation between peak fluorescence response (ΔF/F0) and the number of APs observed in the loose patch recording. The peak fluorescence amplitude was proportional to the number of APs that were counted simultaneously from a loose patch recording (Fig. 1D). This analysis further demonstrated that even a single AP could be detected, as ∼5% fluorescence signal change, well above the background noise of the imaging recording. Because the calcium indicator used in this study has a low dissociation constant (Kd ∼300 nm), calcium ions bound to the indicator molecules took a long time to dissociate; thus, calcium traces did not code multiple APs as multiple sigmoidal responses. Instead, calcium signals showed an integrated, smooth increase in peak value as the number of APs increased. Note that in the following experiments, since we exclusively focus on onset time of the fluorescence change, only the timing of the first spike is taken into account even if there are multiple APs within a burst.
Calcium Imaging and simultaneous extracellular recording. A, Simultaneous loose patch recording (top) and ΔF/F0 (bottom recorded from a spontaneously firing neuron. B, Loose patch and calcium indicator fluorescence traces at higher temporal resolution. Note the close linkage between action potential firing (downward deflections in loose patch recording) and onset of calcium indicator response. C, Plot of calcium indicator response relative to action potential firing (vertical line) derived from simultaneous whole-cell patch and calcium indictor imaging. Note the tight relationship between the indicator response and cell firing. D, Correlation between peak ΔF/F0 value and the number of action potentials in neuronal activity recorded from multiple neurons. (n = 3; Pearson r = 0.956, p < 0.0001).
Spontaneous synchronous network bursts
In most slices prepared from animals aged P4–P6 (15 slices from 11 animals), we observed spontaneous synchronous activity in area CA3. A representative recording from a slice prepared from a P5 rat is shown in Figure 2. Figure 2B shows normalized fluorescence traces (ΔF/F0) for 10 cells, along with a raster plot for 64 of 87 cells that participated in synchronous network activity. Based on the response of these bursts to application of GABA receptor antagonist (described in Effects of blockade of GABAergic inhibition), the frequency of events, and the animal age dependence of occurrence of these events, these network bursts fit the experimental definition of a GDP (Khazipov et al., 1997; Sipilä et al., 2005; Ben-Ari et al., 2007). As is evident in Figure 2, the GDP events occurred repeatedly, i.e., five times in 1 min. There was also occasional, intermittent, nonsynchronous activity independent from the synchronous events in a proportion of neurons. In this particular dataset, 58 of 76 cells (76%) had no activity between synchronous events, and 16 of 76 cells (21%) exhibited isolated firing between synchronous events. The amplitude, or peak fluorescence value, of these nonsynchronous events was often smaller than that seen during the synchronous events. Together with the fact that the ΔF/F0 peak values of most cells during the synchronous events were ∼40% (ΔF/F0), and the fact that the 40% fluorescence response corresponds to multiple number of spikes in Figure 1D, this suggests that neurons fired a burst of APs, i.e., multiple spikes, during synchronous spontaneous network events.
Spontaneous network bursts and onset time analysis. A, An area CA3 steady-state fluorescence (F0) image obtained by fast confocal microscopy of an immature hippocampal slice stained with the calcium indicator OGB-1 AM. B, Synchronous network activity recorded from area CA3. ΔF/F0 traces for 10 cells (top) and raster plot of color-coded fluorescence versus time (bottom) for 64 cells that exhibited synchronous network activity. The synchronous spontaneous network events occurred repeatedly. Note that there is nonsynchronous isolated activity between the synchronous events. C, Comparison of the utility of various image acquisition speeds in capturing neuronal activation onset times during spontaneous network bursts. ΔF/F0 raster plots of grayscale-coded fluorescence change versus time for eight neurons that exhibited spontaneous synchronous activity. Images were captured at 70 ms per frame (C1), a typical speed with full-frame (512 × 512 pixel) resolution, and 11.2 ms per frame (C2), the maximum speed that can be achieved by 4 × 4 binning with our EM-CCD camera. The faster acquisition speed was required to resolve differences in onset times among cells (compare cells with *). D, ΔF/F0 traces (dots) for three different cells within a synchronous network event, together with an onset time estimation by sigmoid curve fitting (line). Dotted vertical lines depict differences in calculated cell onset times for the three neurons. E1, Cellular onset latency distribution of spontaneous network bursts (four bursts). Median onset time within a network burst was defined as zero. If a cell activated earlier than the median value, negative values were assigned to the cell, while positive onset value was assigned to the cells that activated later than median timing. The majority of cells distributed within ∼200 ms, with an SD of 54 ms. E2, Cellular onset latency distribution of evoked synchronous events for the same cells shown in E1. The distribution became narrower, with the majority of cells responding within ∼150 ms, with an SD of 35 ms.
Sufficiently fast acquisition speeds were necessary to measure biologically relevant differences in individual cellular burst onset times during GDP events. To determine the minimal required temporal resolution to capture salient differences in cellular burst onset times, we imaged spontaneous GDPs at two different acquisition speeds; 70 ms per frame, a typical scanning speed for calcium imaging with high spatial resolution (Sasaki et al., 2006, 2008; Allène et al., 2008; Allene and Cossart, 2010) and 11 ms per frame, the maximum speed that could be obtained through 4 × 4 binning of the 512 × 512 pixel EM-CCD camera. Figure 2C shows temporal raster plots of several ROIs from two movies of spontaneous bursts recorded from a slice. The raster plot was zoomed and focused on the rising phase of signal change for a single network burst at both temporal resolutions. Although the two bursts may not be identical, the raster plot in Figure 2C2 resolves differences in cellular burst onset timing to a much greater extent than the recording in Figure 2C1. This demonstrated the requirement for higher temporal resolution recordings to resolve differences in cellular burst latencies underlying spontaneous GDPs.
Onset time distribution
For a series of images taken from the slice, the onset time of synchronous activity was estimated for each ROI (i.e., cell body) by fitting a sigmoidal curve to the rising phase of fluorescence change (see Materials and Methods section). Figure 2D shows an example of onset time estimation for three different cells during a synchronous network event. The sigmoid curves fit the data well. Among 87 cells marked on the image, 15 cells had no synchronous responses and were excluded from the analysis. In addition, several cells had very small responses (<5% ΔF/F0), and the rising phase of these responses was poorly fit (R2 < 0.6). These cells were also excluded from further analysis. After determining the onset time for cells that resulted in a good fit (typically 63–65 cells), as a temporal reference point, the median onset time for the cell population was defined as zero. Cells activating earlier than the median were assigned negative onset time values, while cells that activated later than the median were assigned positive onset times.
Figure 2E1 shows the distribution of onset times for four spontaneous network events recorded from the same slice. As a control experiment, using the same slice, a stimulating electrode was placed in the Schaffer collateral pathway, away from the imaged area, and a burst stimulus (four pulses, 10 μs, 100 Hz, 500 μA) was applied. The electrical stimulus evoked a synchronous network event, and Figure 2E2 shows the onset time distribution of this response. For the spontaneous bursts, onset times ranged from −190 to +286 ms, with an SD of 54 ms. In contrast, for the evoked bursts, the onset time ranged from −207 to +85 ms, with a SD of 35 ms. This suggests that electrical stimulation triggered a different type of network activity, distinct from spontaneous events.
Burst propagation
Since the CA3 layer was tilted slightly to the CCD camera field of view, we rotated the image during analysis to visually align the hippocampal axis, and new image axes were defined as x′ and y′, as shown in Figure 3. Figure 3A shows onset latency as a color-coded two-dimensional map of cell location and cellular burst onset time for a spontaneous network burst. Warmer colors (red) indicate that cells activated earlier than the population median, and cooler colors (blue) indicate that cells activated later. In Figure 3A, red circles (corresponding to cellular ROIs) predominate on the right side of the field, while blue ROIs are located mainly to the left, indicating that cellular activation propagated from right to left, or from CA3b to CA3c, as illustrated in the inset. This can be more clearly visualized by plotting the onset time relative to the cellular x′ location in a scatter graph (Fig. 2, bottom plots). Regression analysis of the cell location versus onset time plot facilitated determination of the direction of burst propagation within the network. In the example in Figure 3A, the fitted latency curve had a minimum on the right side of the image, and the burst took ∼140 ms to travel across the image, thus the propagation speed was estimated to be ∼1.4 mm/s. Note that the burst propagation appeared nonlinear, and that the larger slope in the fitted curve indicates slower traveling speed.
Cellular activation onset analysis of propagation of spontaneous network bursts. A, Top, Color-coded cellular onset time mapping for a spontaneous GDP. Red color indicates that cells are activated earlier than the median timing and blue color indicates cells which activated later. Note that, for this burst, red ROIs are located more to the right side while blue ROIs are located mainly to the left, consistent with propagation from right to left (arrow in lower plot), or from CA3b to CA3c as illustrated in the inset. Note that since the CA3 layer was tilted slightly to the CCD camera field of view, the image was rotated to visually align the hippocampal axis, and new image axes were defined as x′ and y′ as shown in the inset. Bottom, Propagation effect can be seen more clearly by plotting x′ pixel location and onset time in a scatter graph, together with a regression line. B, Top, Onset time plot for a different spontaneous event in the same region. Note that red ROIs are located more to the left side for this burst, consistent with event propagation from left to right (arrow in bottom plot). Bottom, The x′ pixel location/onset time scatter plot for this second burst. Note the opposite propagation direction. See Table 1 for different types of network burst patterns observed in the same slice.
Figure 3B depicts the same analysis method applied to another spontaneous network burst in the same slice. For this burst, red ROIs are concentrated on the left side of the image, while the right side is dominated by blue ROIs, indicating that the burst propagated from left to right, the opposite direction of the burst depicted in Figure 3A. Here the burst traveled in an approximately linear fashion, with a slope of 50 ms per 128 pixels, translating to a speed of ∼4 mm per second. Pearson's correlation analysis between y′ location and onset time resulted in no significant correlation.
In multiple recordings from the slice, we observed a total of 17 spontaneous and 12 evoked network bursts. Based on their characteristics of propagation, we assigned these events to several categories (Table 1). Spontaneous bursts with a slow propagation from right to left (CA3b to CA3c), as depicted in Figure 3A, were seen in the majority of cases (14/17 events) and were defined as a type I burst. A spontaneous burst with faster propagation from left to right (CA3c to CA3b), shown in Figure 3B, was seen only once and was defined as a type II burst. A third burst pattern started in the center of the image and spread laterally in both directions (one event, type III), and last we also observed a disorganized widely distributed onset (one event, type IV). Evoked bursts propagated too quickly to resolve directionality and were classified as type V bursts.
Network burst patterns observed in a slice
Temporal patterns independent of propagation
To begin to resolve additional sources of burst latency variation, we minimized contributions of burst propagation to onset times by restricting our examination to 10 closely neighboring neurons, comparing cellular burst latencies within this CA3 microcompartment. Figure 4A depicts the 10 sampled neurons analyzed, and Figure 4B plots fluorescence traces for two representative network bursts in these cells. Table 2 lists onset time and ranks within these 10 cells for each burst. Based on analysis of burst onset latencies in this anatomically restricted area, cellular activation did not appear to be random. Instead, cells retained a consistent order of activation from burst to burst. Cell 4 repeatedly activated earlier than other cells, while cells 5, 6, and 7 consistently activated later than the majority of cells. To analyze this further, we ranked cells based on relative activation order (shown in Table 2) and the nonparametric ranking order statistic, Kendall's W value, was calculated. Kendall's W is an indicator of ranking order consistency. For these 10 cells and 9 network bursts, Kendall's W was calculated to be 0.63, consistent with a highly repeatable activation order, which was extremely unlikely to occur by chance if activation was based on purely random ordering (p < 0.0001). Based on their activation latencies, we termed cell 4 as an early responder and cells 5, 6, and 7 as followers.
Temporal activation patterns independent of propagation characteristics. Onset time analysis for 10 closely localized neurons, for 9 spontaneous GDPs that exhibited similar propagation patterns (type I in Table 1). A, Fluorescence image with cell id number. B, Fluorescence traces (ΔF/F0) for 10 specified cells for two type I bursts. The traces were ordered, earliest to latest, by the average latency of firing onset during the nine bursts.
Onset Time (and Rank) for 10 closely located cells during GDPs
Based on the above microcompartment analysis, we hypothesized that burst firing order was nonrandom regardless of cellular location (i.e., burst propagation effects), and that individual neurons played distinct roles in defining burst structure (e.g., early responders and followers). To test this hypothesis further, we re-examined our data derived from larger fields, and minimized contributions of burst propagation to assess whether deterministic burst structures were evident in this larger area of the CA3 circuit. We first reanalyzed scatter plots (Fig. 3) for multiple bursts and determined whether individual cells appeared consistently above or below the best fit regression line. If there was a consistent trend that a particular cell fired earlier or later than the timing expected from the propagating characteristic (regression line), this would support the hypothesis that burst structures are deterministic. We minimized the effect of cell location by examining residuals from the spline fit. We then calculated a correlation matrix of these residuals for nine bursts. The median of Pearson's r for all possible 36 pairs was 0.466, and ranged from 0.197 to 0.651 (all but one pair had a significant correlation). This suggested that a consistent cellular activation order was retained, even after location factors were minimized.
Further evidence of deterministic burst structure was observed when comparing the waves that propagated in opposite directions as described in Figure 3, A and B. Pearson's correlation of onset times between two bursts that propagated in opposite directions (Fig. 3A,B) was not significant. However, the correlation of residuals from the spline curves between the two bursts was 0.328 (p = 0.013). This emergence of a correlation between distinct burst types may be additional evidence of a temporal structure of bursts independent of burst propagation contributions.
Latent class analysis
Correlation analysis of residuals from the spline fittings provides support for the existence of a deterministic spontaneous network burst structure. However, correlation by itself does not address the question of whether there are individual cells that are early responders or followers, nor does this type of analysis evaluate the existence of such groups. To begin to identify contributions of individual neurons to burst structure, and establish a statistical framework with which to evaluate the existence of classes (groups) of neurons with differing roles in burst firing order, we applied an LC model, or finite mixture analysis, to explore the clustering of onset time (or onset time residuals) with the specific goal of determining whether there appeared to be a group (class) of neurons that represent early and/or late responders. Details of the LC analysis are discussed in Materials and Methods. Briefly, we fit LC models to the onset times by adjusting the effects of location using a spline. These models allow different numbers of underlying classes (groups) in the data where the membership of each observation to a particular class is unknown and estimated by the procedure. Each model fit to the data differed in the number of classes of cells it assumed. The best fit model or the best number of classes was chosen based on the minimum of an adjBIC.
Figure 5A shows a schematic illustration of the LC model analysis conducted on our dataset. Specifically, onset time was fit to a spline curve in a class-specific manner by only shifting the offset, but not changing the spline function among the classes. By varying the number of classes (from one to seven) and the degrees of freedom for the spline fit (from one to five), adjBIC values were calculated for each scenario, and the scenario expressing the minimum adjBIC value was selected as the optimum model. For our dataset, this procedure identified three classes with a spline of 5 df as the best-fitting latent class model: adjBIC for df = 5: 3187.5 (one class) > 3182.2 (two classes) > 3178.4 (three classes) < 3183.4 (four classes) < 3183.4 (five classes). The three classes were arbitrarily defined as “early,” “middle,” and “late,” referring to the offset of the curve fitting. LC model analysis not only determines the number of latent classes but estimates the posterior probability for each cell belonging to each of the classes. Figure 5B shows the posterior probability of 10 cells as members of three classes shown in bar graphs. For example, cell 1 in Figure 5B had ∼80% chance of being a member of the early class, and ∼20% chance of being a member of the middle class. Cell 3 had a high chance of belonging to the middle class. Cell 4 had chances of being a member of the early (60%) or middle (40%) classes. Thus, class assignment was made based on maximum posterior probability for a cell, with the magnitude of the probability indicating the confidence with which cells were assigned to a class.
Latent class model analysis of cell activation latencies during bursts. A, Schematic illustration of LC analysis with a class-specific linear mixed model. Onset time was fit to class-specific spline curves without changing the shape but shifting intercept in a class-specific manner. Vertical bell-shaped distributions at the right end of each line are class-specific Gaussian distributions. LC analysis was conducted for different number of classes (from one to six); adjBIC values were calculated for each scenario, and the scenario expressing the minimum adjBIC was selected as the optimum number of classes. LC analysis found three classes as the best-fitting model for GDPs. B, Bar graphs of posterior probability of class membership for 10 representative cells (location shown in C) for the three class model, with E (early group), M (middle group), and L (late group). Cells were then assigned to a particular class based on their maximum posterior probability. C, Color-coded class assignment (early, red; middle, green; late, blue) superimposed on the image with the label for 10 cells shown in B. D, Scatter plots of onset time versus x′ location for all cells for nine type I burst responses with the class assignment in color (codes as in C). E, LC analysis of evoked bursts. Note that different groups of cells comprise the early group for evoked bursts compared with spontaneous GDP bursts (C). F, Scatter plots of onset times versus x′ location for all cells for nine evoked bursts with class assignment in color (early, red; major, green).
Table 3 shows a summary of the LC model analysis and class assignments, with the mean response time indicating the difference between the classes. Of 76 cells analyzed, 10 cells (13%) were assigned to the early class, of which 5 cells were assigned with high probability (>90%, 6.5% of total cells). Six cells were assigned to the late class, with the majority of the cells (60 cells) assigned to the middle class. Figure 5C depicts a color-coded class assignment superimposed on the fluorescence image, together with the labeling for 10 cells shown in Figure 5B. Comparing this to Figure 4, the early responders cell 3 and cell 4 in Figure 4 belong to the early group (Fig. 5C, cells 4 and 6), and, similarly, the followers (cells 5, 6, and 7) belong to the late group in LC analysis. Finally, Figure 5D shows scatter plots of onset time for all cells for all nine bursts, with color indicating class assignment. This plot clearly demonstrates that LC analysis is a powerful method for assigning cell classes, considering location factor. The data presented here and above is a representative dataset, containing cells in the early responder group. We have applied a similar LC analysis to spontaneous GDPs from four other slices and found a comparable distribution of early responder cells to that depicted in Figure 5 and Table 3. In total, 5.5 ± 1% of cells (11/199 cells from four animals) were recognized as early responders with high posterior probability (>90%) of class membership by LC analysis.
Class assignment by LC analysis—GDPs and evoked bursts
Temporal patterns of evoked bursts
We also conducted an LC analysis on the network bursts evoked by electrical stimulation (Fig. 2E2; Table 1, type V). The onset time distribution of evoked bursts was narrower than that seen for spontaneous events, as shown in Figure 2E2, and the scatter plots for onset time versus x-pixel location had a best fitting spline (df = 1) with a slope that was not different from 0. This suggested that the propagation speed of these evoked bursts was too fast to be captured and/or the propagation effect was swamped by large variance in the onset times. The correlation between repeated evoked burst measurements (nine bursts) on the same cells from different events ranged from 0.306 and 0.789, with a mean Pearson's r value for all pairs of 0.55 (all significant correlations). Thus, the correlation among repeated measurements collected for the same cell was substantial, which suggests that there is a consistent neuronal firing order, i.e., deterministic burst structure, even in the evoked bursts. The LC model analysis for these evoked bursts (N = 9) generated two classes as the best fit model. The means and SDs of onset times for the two classes are shown in Table 3. We termed the two classes as “early” and “major.” All but seven of the cells were assigned to the major class. Seven cells were assigned to the early class, and, of these, three were assigned with high probability (>90%, Table 3). Figure 5F shows the scatter plot for each cell for nine evoked bursts with class assignment in color. As with the first dataset (spontaneous bursts), the minority class appeared to represent a distinct group based on the difference in the mean values, consistent with the notion of a distinct group of early bursters. Figure 5E depicts a color-coded mapping of class assignment. When compared with the assignment in Figure 5C for spontaneous bursts, it is evident that electrical stimulation triggers network bursts through a different activation pathway, and also that different groups of cells made up the early responder cohort of neurons for these two different bursts, even though they occurred in the same network of cells.
Effects of blockade of GABAergic inhibition
GABAergic neurons play a crucial role in generation of network bursts in immature, early postnatal brains. To investigate the role of GABAergic neurons in generating and shaping the temporal characteristics of spontaneous synchronous bursts, the GABAA receptor antagonist, bicuculline, was applied. There were distinct effects of bicuculline on spontaneous bursts recorded in rats and mice. In rats, application of bicuculline (10 μm) blocked the occurrence of spontaneous synchronous network events as shown in the raster scan in Figure 6A, top (same slice prepared from a P5 rat was recorded in Figs. 2⇑⇑–5). Note that the spontaneous isolated nonsynchronous activity remained intact. The overall occurrence of nonsynchronous activity remained similar before and after applying bicuculline (compare Figs. 2B, 6A, top). Synchronous network bursts could still be evoked by electrical stimulation in the presence of bicuculline, as shown in Figure 6A (bottom), with Figure 6B showing a cumulative distribution of onset times in individual neurons. Interestingly, the onset time distribution for evoked events in bicuculline became wider (Fig. 6B, red) compared with the results obtained for spontaneous GDPs without bicuculline (Fig. 6B, green) (see also Fig. 2E). This suggests that GABAergic inhibition plays a role not only in the initiation, but also in the shaping of synchronous network burst structure. In addition, the number of neurons that were activated by electrical stimulation was reduced in the presence of bicuculline (only ∼50% cells responded to all four electrical stimuli). However, the phenomenon of deterministic activation order was maintained among the responding cells in bicuculline, as the average Pearson's r value for four bursts (six pairs) was 0.613 (range 0.509 to 0.771, all significant correlations).
Effect of blockade of GABAergic inhibition on burst characteristics—two distinct effects. A, Raster plots of color-coded fluorescence intensity (ΔF/F0) versus time for 50 cells, in the presence of bicuculline (top) and during stimulation in the presence of bicuculline (bottom, 15 s intervals). Addition of bicuculline (10 μm) completely blocked GDP bursts, recorded from the same slice as shown in Figure 1⇑–3, prepared from P5 rat. Note that bicuculline did not block nonsynchronous activity and that synchronous bursts could still be induced by electrical stimulation. B, Plot of onset time distribution for GDPs (green), evoked responses in the absence of bicuculline (blue), and evoked responses in the presence of bicuculline (red). Note that evoked responses without bicuculline exhibited a tighter temporal distribution and bursts evoked in the presence of bicuculline a wider temporal distribution than GDPs. C, Bicuculline did not block spontaneous network bursts in every preparation. Top, Raster plots of color-coded ΔF/F0 versus time for 75 cells for spontaneous synchronous network bursts (P5 mouse). Bottom, 10 μm bicuculline did not block spontaneous network bursts but induced extremely synchronous and intense events. These network bursts were termed interictal bursts in the presence of bicuculline. D, Plot of cellular activity onset time distribution of interictal bursts (purple dots) revealed that cellular activation times were much tighter (30–40 ms window) than spontaneous bursts without bicuculline (green; 200–250 ms window). E, Fluorescence traces (ΔF/F0) from three representative cells in C. Spontaneous bursts (GDPs) on the left and spontaneous bursts under bicuculline (interictal bursts) on the right. Note that right traces in each panel are at an enhanced time scale. The intense interictal activity resulted in higher peak ΔF/F0 values in individual neurons, which decayed more slowly than that seen during GDPs. The rising phase of ΔF/F0 traces during interictal activity was also faster than that of GDPs.
Interictal bursts
We observed synchronous network bursts that were completely blocked by bicuculline in many slices, particularly those prepared from rats. However, bicuculline did not always block spontaneous synchronous bursts. In 11 slices prepared from 7 mice, aged between P4 and P6, bicuculline application decreased the frequency of occurrence of spontaneous synchronous events, and increased the amplitude and intensity of these responses. This observation was also described by Le Magueresse et al. (2006), who attributed this phenomenon to the emergence of epileptiform interictal discharges upon GABA blockade. Figure 6C illustrates this effect of bicuculline, depicting time raster plots of cellular fluorescence responses during spontaneous network bursts before (Fig. 6C, top) and after (Fig. 6C, bottom) the addition of 10 μm bicuculline to a slice prepared from a P5 mouse. Interestingly, and in contrast to the evoked burst responses in bicuculline depicted in Fig. 6A, the onset time distribution of bicuculline-transformed interictal activity was much tighter than that seen for spontaneous events before treatment (Fig. 6D), with all activation latencies falling within a range of three frames (33 ms), compared with >200 ms for control slices. The intense interictal activity resulted in higher peak ΔF/F0 values in individual neurons (peak ΔF/F0 up to 80%), which decayed more slowly (∼5 s) than that seen during GDPs (which exhibited a peak ΔF/F0 up to 40% with a decay in ∼2 s). The rising phase of ΔF/F0 traces during interictal activity was also faster than that of GDPs (Fig. 6E).
To better resolve differences in cellular activation latency during this very tightly distributed interictal activity, we enhanced our imaging speed through the use of a faster CCD camera, which allowed us to increase our image acquisition rate to 350 Hz. The onset time for interictal events varied within a range of 45 ms (16 frames at enhanced resolution) as shown in Figure 7A. About one-third of cells had spontaneous nonsynchronous runs of responses (data not shown), which on occasion could occur immediately before a synchronous interictal burst, skewing the onset time analysis. In these datasets from six interictal bursts, four responses from four cells had responses up to 100 ms earlier than the majority of responses. Because this was not a repeatable occurrence in individual cells, and spontaneous responses could also occur not in association with network bursts, these four data points were discarded as outliers. This left responses from a total of 37 cells that had responses during all six waves. Correlations between repeated measures on the same cell from different bursts ranged from −0.26 to 0.55, with a median correlation of 0.06 (p value not significant). Thus, the correlations were substantially weaker than for the first dataset (Fig. 5). This indicates that interictal bursts exhibited a less deterministic burst structure than that seen during GDPs, associated with a more random neuronal activation order in repeated bursts.
Interictal bursts. A, Onset time distribution of interictal network bursts (imaged at 350 Hz). The cell activity latency clustered tightly (within ∼40 ms), with four outlier values (arrows). B, Latent class model analysis resulted in two classes as the best fit model (see Table 4). Median onset time for each cell is plotted together with the range of values for six interictal bursts, with class assignment in color. The four most extreme negative outcomes in the dataset were outliers, and are evident in the wide range for one majority (green) and three minority (red) class cells. The median cell latency values for all cells overlapped.
To further examine this issue, we applied the latent class model to interictal burst discharges, using data from all 49 cells with a response to one or more of the six bursts, including the outliers. The latent class analysis identified two classes as the best fit model. Means and SDs for the two classes from the best-fitting latent class model are shown in Table 4. We termed the classes as “minor” and “major”. All but three of the cells were assigned to the major class. Three cells were assigned to the minor class. Of these, two cells were assigned to the minor class with high posterior probability (>90%). In contrast to the data in Figure 5 derived for GDPs, the three cells in the minor class did not represent a distinct group, as the median onset time for these cells was not different from that in cells in the major class (Fig. 7B). Rather, each cell in the minor group had three of the four most extreme negative outcomes in the multiburst dataset, i.e., outliers. These outliers appear to be highly influential in forcing the latent class model to have two groups. Thus, while our analysis found two groups in these data, the results are much less consistent with the notion of a group of early responders. Our conclusion from this analysis is that interictal bursts exhibit little in the way of robust class structure, and therefore are composed of a population of neurons exhibiting stochastic (nonordered) firing latencies during synchronous network activation.
Class assignment by LC analysis—interictal bursts
Discussion
GDP network bursts in CA3 are deterministic
In the present study, derived from a detailed microcircuit analysis, we found that there is a consistent, deterministic neuronal activation order underlying GDPs, which is distinct from the effects of cellular location or the propagation characteristics of these events. Statistical approaches demonstrating this phenomenon included a high correlation coefficient in neuronal firing latency from burst to burst after minimizing the cellular location factor, and distinct classes representing both early and late firing neuron populations in LC analysis. This contrasted with cell firing orders resolved during spontaneous interictal bursts, which appeared stochastic, with no cells repeatedly leading or lagging in firing latency during these events.
The spontaneous synchronous network bursts shown in Figures 2⇑⇑–5 fit the characteristics of GDPs reported in the literature: they are observed in early postnatal ages, recur at repeated intervals, are sensitive to GABAA receptor antagonists (bicuculline), and exhibit a broad temporal spread of neuronal activation timing (Ben-Ari et al., 2007). GABA operates primarily via chloride-permeable GABAA receptor channels. During early development, neurons have a higher intracellular chloride concentration due to immature expression of transmembrane chloride transporters. Because of this, GABA triggers an efflux of chloride, membrane depolarization, summing in a net excitatory effect promoting AP firing. GDPs are generated by synergetic interactions between recurrent excitatory GABAergic and glutamatergic synapses, combined with intrinsic regenerative conductances of CA3 neurons (Khazipov et al., 1997; Sipilä et al., 2005; Ben-Ari et al., 2007). Factors contributing to delayed action potential firing in CA3 neurons during GDPs have been subject to detailed analysis (Valeeva et al., 2010), where it was suggested that neuronal responses during GDPs were initiated by depolarizing IPSPs driven by GABAergic interneurons. These IPSP depolarizations were subthreshold for action potential firing. This required activation of persistent Na channels by the IPSP to further depolarize the neuron, which initiated delayed AP firing.
In our study, we found that neurons fired APs in a distribution spanning 200–300 ms during spontaneous GDP network bursts, even though, in many cases, these cells were close neighbors. Given that interneuron collaterals frequently distribute in intensely local arbors, this appears inconsistent with the contention that GABAergic interneuron activity alone can control the precise temporal activation order of multiple, adjacent pyramidal cells. Rather, we hypothesize that the activity of GABAergic neurons may control the global propagation direction, creating conditions amenable for pyramidal cells to fire APs, but the activation timing is determined by the distinct intrinsic properties, and perhaps anatomy, of the pyramidal cells (Ropireddy et al., 2011). We further hypothesize that the early, middle, or late class membership determined by the LC analysis is directly related to intrinsic, cellular properties such as transmembrane chloride gradients, resting membrane potential, and/or sodium or calcium channel activity. This may allow subsets of cells to respond more strongly (activate more readily) to smaller depolarizing inputs. We also speculate that class membership determined by the LC analysis could also be related to differences in synaptic connectivity or dendritic anatomy of individual pyramidal neurons. Neurons with more widely distributed dendrites may collect synaptic input (activity) from broader areas relative to a lesser neighbor, which could allow this cell to become a watchful early responder in wave propagation. Similarly, if a neuron has a more tightly distributed dendritic tree, it may be able to collect recurrent collaterals only from its own immediate neighbors, relegating this neuron to a follower role.
The fact that waves propagated in different directions within the same microcircuit (Fig. 3) suggests that there are numerous burst-initiating locations within a slice. As depicted in Fig. 5D, our analysis showed that early responders were not necessarily the earliest activating cells within an image field, suggesting these early responders were not burst initiators. Because of their location (cell body layer) and cell body morphology (no larger than other cells), and the fact that we did not observe continuous firing in these early responders (which might be a characteristic of interneurons), we speculate these neurons are subsets of CA3 pyramidal cells. Early responders, once activated, could activate other principal neurons or interneurons through recurrent collaterals. Those secondary neurons could further activate other principal cells. Thus, it is possible that, without those early responders, GDPs may not propagate as effectively, nor recruit the same large proportion of neurons into the burst.
If the role of GDPs is to form functional connectivity by fire together-wire together mechanisms in developing hippocampal networks (Ben-Ari et al., 2007), the deterministic firing order we saw in our experiments (Fig. 4) may have a significant role in forming temporally important wiring patterns. The fact that GDPs are so synchronous and activate such a large proportion of the network makes it unlikely that the GDP itself can carry cognitive information. However, the temporal patterns activated during the GDP may be important in future coding, where any two neurons that express associated temporal delays can carry important information in a temporal coding mechanism (VanRullen et al., 2005). It is also possible that the delayed activation in CA3 may have significant effects on downstream (CA1) signal transmission by delayed coactivation mechanisms (Izhikevich and Hoppensteadt, 2009).
The importance of deterministic firing order in synchronous network bursts and the possible role of this phenomenon in neuronal coding has been examined previously (Beggs and Plenz, 2004). These investigators studied cultured cortical slices using a multielectrode array and analyzed network bursts. They also found several repeatable activation orders, albeit within a time scale of ∼10 ms (compared with the 200–300 ms timescale of GDPs) (Fig. 4). In another multielectrode array recording study from sensory cortex in vivo, population spikes caused by sensory stimuli as well as spontaneous events were analyzed (Luczak et al., 2009). These investigators also found a repeatable, sequential neuronal activation order that lasted ∼100 ms. They argued that the patterns of activation order and frequency response comprised the “vocabulary” of coding for sensory inputs.
Bicuculline-induced interictal bursts are stochastic
Without the inhibitory synaptic transmission mediated by GABAergic neurons and possibly GABA release from mossy fibers (Safiulina et al., 2010), the CA3 network is mainly composed of glutamatergic excitatory principal cells and their recurrent collaterals. In contrast to GDPs, our experimental analyses could not identify any evidence for ordered burst structure (e.g., early or late responder groups) in interictal bursts occurring in the presence of bicuculline, seen exclusively in slices prepared from mice. Interictal bursts appear to be quite distinct from GDPs, implicating inhibitory synaptic transmission as critical both in the initiation and shaping of GDPs. Several studies have examined network burst generation mechanisms in disinhibited hippocampal networks (Miles and Wong, 1983; de la Prida et al., 2006; Jones et al., 2007). Using electrophysiology, de la Prida et al. (2006) demonstrated that initiation of spontaneous population bursts in disinhibited CA3 region followed a buildup of synaptic excitation and firing to a threshold level, which could take as long as 50 ms. This group also showed that activation of a single neuron could induce a network burst when the timing matched with the recovery period from the previous bursts (Jones et al., 2007) (see also Miles and Wong, 1983), and furthermore, that ∼50% of cells they investigated had the capability to induce a network burst, suggesting that many cells could initiate buildup of cellular activity that leads to a population burst, supporting a stochastic process underlying burst generation. In our experiments, it is possible that the outliers we saw in Figure 7B might correspond to a similar buildup of activity. Since the cells that fired early did not retain this role repeatedly from burst to burst, we viewed the early firing phenomenon to be a stochastic process. The median onset time of these cells was not different from the remaining population (Fig. 7B), suggesting that the network was composed of a relatively homogeneous population. This may be due to the fact that blocking interneuronal contributions to spontaneous network bursts created a much simpler network composed of only glutamatergic recurrent circuits.
Computational modeling studies examining mechanisms of network burst generation in homogeneous networks composed of only one cell type with varied connectivity and excitability parameters have found that a homogeneous population could generate network bursts when cells have high connectivity and random wiring (Izhikevich, 2003). Additional modeling of neuronal avalanches showed that synchronous bursts were generated in noisy systems and did not require a self-organized, deterministic network (Benayoun et al., 2010). These modeling studies suggest that network bursts such as interictal activity in our experiments could be generated by stochastic mechanism without having a specific cell type that initiates these events.
In contrast to the above discussion on homogeneity, populations of bursters have been found in area CA1 in the high K+ model of epilepsy as well as in tissue prepared from animal models of epilepsy (Jensen and Yaari, 1997; Sanabria et al., 2001; Becker et al., 2008; Chen et al., 2011). These investigators hypothesized that CA1 population bursts were generated by excitatory input from CA3 activating a burster population in CA1 first, followed by these bursters recruiting regular pyramidal cells by local electrical field or local synaptic interactions. In area CA1 of brain slices prepared from epileptic animals, a significant proportion of cells (up to 47%) were found to be bursters, and activation of an individual burster neuron could trigger interictal events. The computational modeling study by Bogaard et al. (2009) also explored the role of heterogeneous populations on network burst generation. When two types of bursting characteristics were mixed, only a small percentage (10%) of a bursting-favorable cell types could have a significant effect on network burst generation and its pattern. This complex burst-generating mechanism, where heterogeneous populations (deterministic, defined by intrinsic properties of individual neurons) act on a stochastic process—a hybrid mechanism—may be the realistic model of actual bursting mechanism in epileptic brains.
Notes
Supplemental material for this article is available at https://dbe.med.upenn.edu/biostat-research/MaryEPutt. Latent class analysis. Thi material has not been peer reviewed.
Footnotes
This work was supported by grants from the NINDS.
- Correspondence should be addressed to Dr. Douglas A. Coulter, The Children's Hospital of Philadelphia, 3615 Civic Center Boulevard, Abramson Research Center, Room 410D, Philadelphia, PA 19104-4318. coulterd{at}email.chop.edu
