Abstract
Ongoing interactions among cortical neurons often manifest as network-level synchrony. Understanding the spatiotemporal dynamics of such spontaneous synchrony is important because it may (1) influence network response to input, (2) shape activity-dependent microcircuit structure, and (3) reveal fundamental network properties, such as an imbalance of excitation (E) and inhibition (I). Here we delineate the spatiotemporal character of spontaneous synchrony in rat cortex slice cultures and a computational model over a range of different E–I conditions including disfacilitated (antagonized AMPA, NMDA receptors), unperturbed, and disinhibited (antagonized GABAA receptors). Local field potential was recorded with multielectrode arrays during spontaneous burst activity. Synchrony among neuronal groups was quantified based on phase-locking among recording sites. As network excitability was increased from low to high, we discovered three phenomena at an intermediate excitability level: (1) onset of synchrony, (2) maximized variability of synchrony, and (3) neuronal avalanches. Our computational model predicted that these three features occur when the network operates near a unique balanced E–I condition called “criticality.” These results were invariant to changes in the measurement spatial extent, spatial resolution, and frequency bands. Our findings indicate that moderate average synchrony, which is required to avoid pathology, occurs over a limited range of E–I conditions and emerges together with maximally variable synchrony. If variable synchrony is detrimental to cortical function, this is a cost paid for moderate average synchrony. However, if variable synchrony is beneficial, then by operating near criticality the cortex may doubly benefit from moderate mean and maximized variability of synchrony.
Introduction
Neurons whose dynamics are synchronized are thought to play an important role in brain function, because their combined influence on other neurons is greater than that of asynchronous neurons (Kelso, 1995; Varela et al., 2001; Kopell et al., 2010; von der Malsburg et al., 2010). Abnormal synchrony measured during ongoing, spontaneous activity can indicate dysfunction of the underlying cortical network. For example, excessive synchrony occurs during epileptic seizures (Steriade, 2003; Garcia Dominguez et al., 2005) and Parkinson's disease (Boraud et al., 2005; Uhlhaas and Singer, 2006), while abnormally weak synchrony has been associated with disorders such as schizophrenia (Spencer et al., 2003; Uhlhaas and Singer, 2010) and autism (Wilson et al., 2007; Yizhar et al., 2011). Thus, a healthy cortical network must be regulated such that spontaneous synchrony avoids pathological extremes and remains moderate on average.
In addition to its average level, the variability of spontaneous synchrony is also likely to influence cortex function. Transient periods of synchrony between two groups of neurons can serve as a temporal “windows of opportunity” during which interactions between them are more effective (Varela et al., 2001; Womelsdorf et al., 2007; Canolty et al., 2010). In this view, variability of ongoing synchrony fluctuations entails dynamically changing windows of opportunity. Such variability may explain trial-to-trial variability of response to input from external sources (Arieli et al., 1996; Poulet and Petersen, 2008; Cafaro and Rieke, 2010). In the absence of external stimuli, spontaneous activity occurs in spatiotemporal patterns which often reflect previous experience (Ohl et al., 2001; Kenet et al., 2003; Ji and Wilson, 2007; Han et al., 2008) and, thus, may influence consolidation and maintenance of memory. Therefore, a large and variable repertoire of spontaneous synchrony may be linked with a corresponding repertoire of experiences (Luczak et al., 2009; Berkes et al., 2011).
Both the mean and variability of spontaneous synchrony are expected to depend on intrinsic cortex properties, like the relative influence of excitatory (E) and inhibitory (I) neuronal populations, which is the focus of our study. In disorders such as epilepsy, an imbalance favoring excitation may underlie the excessive synchrony of a seizure (Dichter and Ayala, 1987). However, a more general understanding of how cortical E–I affects synchrony dynamics remains lacking. Here we studied spontaneous phase synchrony in rat cortex slice-cultures and a computational model over a range of different E–I conditions. Our computational model predicted that moderate mean and high variability of synchrony would emerge together for a particular balanced combination of E and I called “criticality.” Criticality can have different definitions in different contexts. Here we refer to the definition used in the study of non-equilibrium phase transitions (Hinrichsen, 2006). At criticality, the network dynamics take the form of neuronal avalanches (Plenz and Thiagarajan, 2007; Shew et al., 2009). Similar phenomena were previously predicted for coupled oscillators operating near criticality (Daido, 1990). Our experiments confirmed these predictions; moderate mean and peak variability of synchrony were found uniquely under the E–I condition at which we observed neuronal avalanches.
Materials and Methods
Organotypic cortex cultures.
Coronal slices from rat somatosensory cortex (350 μm thick) and the midbrain (VTA; 500 μm thick) were taken from newborn rats (postnatal day 0–2; Sprague Dawley, of either sex) and cultured on a poly-d-lysine-coated 8 × 8 multielectrode array (MEA; MultiChannel Systems; 30 μm electrode diameter; 200 μm interelectrode distance, 59 electrodes used for recording). This coculture system captures many aspects of the development of superficial cortical layers observed in vivo (Gireesh and Plenz, 2008) and has been described in detail previously (Gireesh and Plenz, 2008; Shew et al., 2009, 2011). In short, the slice cultures were grown inside a sterile, closeable chamber attached to the MEA. Individual MEAs were attached to a recording head stage inside an incubator (MEA1060 with blanking circuit; ×1200 gain; bandwidth 1–3000 Hz; 12 bit A/D in range 0–4096 mV; MultiChannel Systems). The recording setup allowed for repeated measurements under stable, sterile conditions from single cultures for weeks. Spontaneous local field potential (LFP; 4 kHz sampling rate; measured against a common reference electrode inside the bath) was obtained from 1 h of continuous recordings of extracellular activity, and low-pass filtered off-line with a cutoff at 50 Hz (phase neutral, fourth-order Butterworth).
Pharmacology.
Bath-application of antagonists for fast glutamatergic or GABAergic synaptic transmission was used to change the excitability of the cortical network. The normal condition (no-drug) followed by a drug condition for each culture was typically studied within a short time (∼3 h) to minimize potential non-stationarities during development. Stock solutions were prepared for the GABAA receptor antagonist picrotoxin (PTX), the NMDA glutamate receptor antagonist AP5, and the AMPA glutamate receptor antagonist DNQX. Working solutions were obtained by adding 6 μl of these stock solutions to 600 μl of culture medium in the MEA chamber to reach the following final drug concentrations (in μm): 5 PTX, 20 AP5, 10 AP5 + 0.5 DNQX, 20 AP5 + 1 DNQX. We intentionally chose drug concentrations which perturb synaptic transmission without causing total blockade of receptors. For example, at a working concentration of 5 μm, PTX is known to reduce the amplitude of IPSCs by 30–75%, while total blockade requires ∼50 μm (Sedelnikova et al., 2006; Korshoej et al., 2010). Furthermore, we previously showed that carbonoxolone has negligible effects on cortical dynamics in our in vitro model for the developmental time period considered here (Gireesh and Plenz, 2008). Thus, we expect the dynamical changes observed not to be confounded by a relative change in the contribution of gap-junction versus synaptic transmission related coupling. After measurement under each drug condition, the culture was washed by replacing the culture medium with 300 μl of conditioned medium mixed with 300 μl of fresh, unconditioned medium. Conditioned medium was collected from the same culture the day before drug application. Most cultures recovered to their predrug condition within ∼24 h following wash.
Bursts.
First, we identified negative peaks in the LFP that fell below −4 SDs of the electrode noise. These point events were called nLFPs. Electrode noise is defined for each electrode during quiescent periods with no bursts. We previously showed that nLFPs are correlated with increased spiking in the cortical cultures (Gireesh and Plenz, 2008; Shew et al., 2009). For each nLFP, we obtained three quantities: timestamp, electrode on which it occurred, and amplitude. We then defined bursts as spatiotemporal clusters of nLFPs within which consecutive nLFPs (from any electrode) were separated by less than a time τ. Bursts were typically comprised of many nLFPs covering many electrodes and lasting from a few milliseconds to hundreds of milliseconds. Following previous studies (Shew et al., 2009, 2011), the threshold τ was chosen to be greater than the short timescale of inter-nLFP intervals within a burst, but less than the longer timescale of interburst quiescent periods (τ = 86 ± 71 ms for all cultures). For the purposes of computing κ (see below), the size s of a burst was quantified as the absolute sum of all nLFP amplitudes within the burst. Thus, the size s captures three aspects of the burst: spatial extent, duration, and the amplitude of the nLFPs. We also studied the first two of these properties separately. The spatial extent of a burst, a, was defined as the number of electrodes from which an nLFP was recorded during the burst. The duration of a burst, d, was defined as the interval between the first and last nLFP within the burst.
Definition of neuronal avalanches and κ.
Here and previously, neuronal avalanches are defined by a probability density function (PDF) of burst sizes that is a power-law with exponent close to −1.5, i.e., Pr(s)∼s−1.5 (Beggs and Plenz, 2003; Shew et al., 2009; Klaus et al., 2011). This definition was motivated by theory of critical phenomena (Hinrichsen, 2006). One misconception is that any large burst of neural activity is a neuronal avalanche. On the contrary, one must measure many bursts, compute the burst size PDF, and finally, only if the PDF is of the form Pr(s)∼s−1.5 can one conclude that the measured bursts were neuronal avalanches.
Under different experimental conditions, one may encounter forms of burst size PDFs that differ from the PDF expected for neuronal avalanches. To compare and quantify how far the observed burst activity deviates from avalanche dynamics, we use a previously defined measure κ, which assesses how close an experimentally measured burst size PDF is to Pr(s)∼s−1.5 (Shew et al., 2009). If the size is bounded between a minimum size l and a maximum size L, then the properly normalized distribution is Pr(s) = Cs−1.5, where the constant C = 0.5
where the βk values are m = 10 discrete burst sizes logarithmically spaced between the minimum and maximum burst size observed in the experiments. Thus, if an observed burst size PDF is an exact −1.5 power-law, then κ = 1. Previous work with computational models has shown that κ = 1 can occur at the critical point of continuous phase transition, at which neuronal avalanches are expected to occur, consistent with the universality class of directed percolation (Shew et al., 2009, 2011). In this context, κ > 1 indicates supercriticality and κ < 1 indicates subcriticality.
Our measurements of burst area, duration, synchrony, and entropy (synchrony and entropy are defined below) were plotted against κ for all n = 47 experiments. Experiments with similar κ were averaged together by binning κ linearly in steps of 0.15 from 0.55 to 1.75.
Phase synchrony.
We first obtained a phase trace from each LFP trace. This approach treats the LFP signal f(t) as the real part of a complex signal, called the analytic signal z(t). The imaginary part of z(t) is the Hilbert transform of f(t), and z(t) = f(t) + iH[f(t)]. We implemented this approach using the MATLAB (MathWorks) function hilbert. The phase trace is calculated as:
In this way, we extracted the instantaneous phase of every electrode over the entire duration of the recording. For pure periodic signals, e.g., f(t) = cos(ωt + θ0), this approach yields the phase θ(t) = ωt + θ0, where θ0 is the phase at t = 0. The phase is insensitive to amplitude, e.g., the phase traces are identical for Acos(ωt) and Bcos(ωt), where A ≠ B. Note that in our analysis we “wrapped” the phase so that it was always between −π and π. Similarly, for aperiodic signals like our LFP measurements, the instantaneous phase θ(t) is related to zero crossings and peaks/troughs of the LFP signal, while the time derivative dθ(t)/dt gives an estimate of the instantaneous frequency of the LFP.
Phase and phase difference histogram.
To display the dynamics of phase and phase synchrony in a way that is amenable to visual examination, we used two analysis techniques. First, the dynamic phase histogram displays phase dynamics for all 59 electrodes simultaneously. At each time point, a histogram of the 59 phases from −π to π (resolution 0.05π) is computed and displayed as one vertical strip. A color code indicates the number of electrodes in a given phase bin. For example, red pixels indicate the times when many electrodes have similar phase, i.e., are phase synchronized. A second method to visualize phase synchrony is the dynamic phase difference histogram. We first compute the difference in phase for every pair of electrodes (n = 59 · (59 − 1)/2 = 1711 pairs). A histogram of the phase differences is computed and displayed as one vertical strip at every time point (resolution 0.05π) and color coded. Phase synchrony is apparent as times when many pairs have near-zero phase difference. We note that some previous studies use the term “phase synchrony” more generally to describe any constant phase difference, not just zero phase difference. This visualization approach is well suited to identifying such dynamics. In our experiments, we did not observe non-zero phase difference synchrony with the exception of a tendency for anti-phase locking, i.e., near-π phase differences, for high κ, i.e., under disinhibited conditions, 1.2 < κ < 1.6. Hereafter we shall refer to near-zero phase difference as phase synchrony and near-π phase difference as anti-phase locking.
Anti-phase locking.
To quantify the degree of anti-phase locking, we selected 8 experiments for each of the three conditions κ < 1 (disfacilitated), κ ∼ 1 (normal), and κ > 1 (disinhibited). From each experiment, we randomly chose 100 bursts. For each burst, we computed the time-average of the dynamic phase difference histogram. In this calculation, we included a 5 ms window preceding and following each burst, which is necessary to capture residual synchrony around burst occurrence. Next, we averaged histograms across bursts and computed the difference between the count of electrode pairs with phase difference near π, i.e., within π ± 0.05π, and the minimum of the histogram. This value which quantifies the degree in anti-phase locking is reported in the variable APL in Results. In addition, we studied whether anti-phase locking tended to occur between electrode pairs that spanned superficial and deep cortical layers. The midline of the MEA approximately divides superficial layers from deeper layers in the cortex cultures (Gireesh and Plenz, 2008) and we computed the fraction of anti-phase-locked pairs which spanned the midline. For randomly evolving phases, this number would be 0.508, which is the fraction of all possible electrode pairs that span the MEA midline.
Network synchrony and burst synchrony.
A simple way to quantify the network-level phase synchrony as a function of time is provided by the Kuramoto order parameter, r(t), defined as
where n = 59 is the number of recording electrodes among which synchrony is estimated (Strogatz, 2001; Arenas et al., 2008). We were particularly interested in the synchrony during burst activity, which we quantified using the following three measures. First, the “network synchrony” of the kth burst, SNk, was defined as
where tk is the burst start time and dk is the burst duration. Second, we calculated the “instantaneous network synchrony” for the kth burst, SINk, by normalizing SNk by burst duration dk,
We note that SINk is now bounded between 0 and 1. Third, we computed the “instantaneous burst synchrony,” SIBk, by applying Equations 3–5 only to those electrodes that were active during the kth burst. More precisely, r(t) was replaced by rm(t) computed only among the set of m sites, Ek, which were active, i.e., displayed an nLFP, during the kth burst
Note that rm(t) represents a measure of within burst synchrony normalized by the area of the burst, i.e., number of electrodes that participate. RC(m) is a correction to account for the expected level of noise for m sites with randomly evolving phases. RC(m) was estimated by direct numerical simulation (104 time steps) of rm(t) for m random phases. By subtracting RC(m), the actual values of rm(t) are comparable even if m is different for different bursts. This correction was not applied to SINk or SNk and if it was applied it would not change any conclusions because these quantities were always computed with the same number of sites. SIBk is bounded between −1 and 1 because both terms in Equation 6 are bounded between 0 and 1.
The averages S̄N, S̄IN, and S̄IB were obtained by averaging over bursts. For example,
where M is the total number of bursts recorded in an experiment.
Entropy measurements.
In our results, we present entropy measurements of burst spatial extent, burst duration, and the three types of synchrony defined above. First, we computed the quantity in question, say x (which can be burst area, burst duration, SNk, SINk, or SIBk), for all M measured bursts. Then, a probability distribution of x was estimated; the probability pi that x fell within the range bi < x ≤ bi+1 was estimated as the number of x values in that range divided by M. The entropy of the x distribution was computed as
where B is the number of bins used to make the distribution. Since entropy is sensitive to the choice of bins, the bin divisions were fixed for comparisons across experiments, i.e., for different κ. Importantly, we found that the peak in entropy persisted for a very wide range of numbers of bins, from B = 4 to 4000 bins. However, the numerical values of entropy in terms of bits are not meaningful, except perhaps in reference to log2(B), which is the maximum possible entropy computed from any distribution with B bins.
Burst area is a discrete variable ranging from 1 to 59, thus, we used one bin for each possible burst area. In contrast, at a temporal resolution of 0.25 ms, burst duration was treated as a continuous variable and we constructed distributions using 100 logarithmically spaced bins ranging from the shortest to the longest duration observed from all experiments (∼10 ms to several seconds). Likewise, the distribution used to compute the entropy of SN was created with 100 logarithmically spaced bins spanning the range of observed SN values. For the bounded quantities, SIN and SIB, 100 linearly spaced bins between 0 and 1, −1 and 1 were used respectively. We note that negative values can result for SIB due to the RC(m) correction.
Finally, entropy can be biased toward low values if a small number of samples are used to estimate the probability distribution. To account for this bias, we estimated subsampling corrections to entropy following an established strategy (Magri et al., 2009) and found that in all cases, the entropy corrections were small compared with the variation from experiment to experiment. Thus, subsampling bias does not significantly impact our conclusions.
Power spectra.
Power spectra were computed via Welch's method using the MATLAB function pwelch. First a power spectrum was computed for each burst. A period of 250 ms preceding and following each burst was included in the calculation and only bursts of duration >10 ms were included. Thus, the spectrum of each burst was computed using time traces at least 510 ms in duration, but typically longer. Note that we did not split the duration of the burst into 8 segments as would be the default option of MATLAB's pwelch function. Instead, we split it into 375 ms windows with 187.5 ms overlap. The FFT was computed on each 375 ms segment and then averaged over segments. Thus, the minimum frequency in our spectra and the resolution was ∼3 Hz. Finally, all spectra from all bursts were averaged and normalized by maximal power to facilitate comparison across experimental conditions.
Statistical analysis.
For determining the statistical significance of differences in anti-phase locking for different drug conditions, we first used a one-way ANOVA to establish that at least one κ category was different from at least one other. Next we performed a post hoc test of significant pairwise differences between the κ categories using a t test with the Bonferroni correction for multiple comparisons.
Computational model.
Our computational model consisted of 59 groups of neurons, which were intended to represent the 59 recording sites in our experiment. The number of neurons per group/site was varied from 30 to 80 (in steps of 10), resulting in network sizes ranging from Q = 1770 to 4720 total neurons. Every group contained 20% inhibitory and 80% excitatory neurons. A Q × Q matrix specifying the connection strength wij between each pair of neurons was constructed as follows. First, all connections were drawn randomly from a uniform distribution on [0, 1]. Next, 20% of the matrix columns were multiplied by −1; these are the outgoing connections from the inhibitory neurons. Then, the entire matrix was multiplied by a constant such that the average matrix entry was exactly equal to 1/Q. This matrix was defined as the “normal condition” and the entries of the matrix were perturbed away from this baseline to model disinhibited and disfacilitated conditions as detailed below. It is a network that is all-to-all connected with inhomogeneous, directed connections with mean and SD of order Q−1. The dynamics of the network was updated synchronously as follows,
where the state bi of the ith neuron was either 0 (quiescent) or 1 (spiking), Θ[·] is the Heaviside step function, wij is the entry of the connection matrix specifying the connection from neuron j to neuron i, and ζ is a random number from a uniform distribution on [0,1] drawn afresh for each neuron update. In short, these are binary, probabilistic, integrate-and-fire neurons. A refractory period during which no spikes could occur was enforced (overruling Eq. 9) for 4 time steps following each spike.
Bursts of network activity were modeled one at a time, always initiated as a single randomly chosen active excitatory neuron and allowed to evolve until either the activity ceased or until 100 timesteps were reached. A total of 1000 bursts were modeled for each condition. The size of a burst was defined as the total number of spikes that occurred during the burst. We note that in previous work we showed that spike count during a burst was proportional to the burst size definition based on sum of nLFP (Shew et al., 2009). Thus, defining burst size in terms of spike count is appropriate for the computational model. Size distributions of these bursts were used to parameterize the model dynamics in terms of κ as defined above. Moreover, each condition was modeled on 60 different connection matrices (10 different networks of each size) to represent the possible variability from one culture to another in the experiments. The error bars in the results represent variation among the dynamics of these 60 different networks.
To model the pharmacological manipulations made in our experiments, we simply multiplied either the positive or the negative entries of the connection matrix by a constant between zero and 1, modeling the effects of either disfaciliation (AP5/DNQX) or disinhibition (PTX) respectively. Since the average wij value of the baseline connection matrix was 1/Q, disfacilitation or disinhibition resulted in a mean wij <1/Q or >1/Q, respectively. As average wij is increased from <1/Q to >1/Q, the dynamics of this computational model undergo an abrupt change at the critical point specified by 1/Q. At the critical point the computational model generates neuronal avalanches, i.e., the burst size distribution is a power law with exponent near −1.5. Our computational model is similar to others which have been used to study criticality in neural networks (Beggs and Plenz, 2003; Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006; Larremore et al., 2011), but is different from these because we include inhibitory neurons.
Before analyzing the computational model data, aggregate population signals PN(t) for each model site were computed PN(t) =
Results
Neuronal burst area and duration have moderate mean and maximum entropy near κ = 1
We studied ongoing network activity recorded in organotypic tissue cultures grown on planar integrated MEAs. Following established techniques (Gireesh and Plenz, 2008; Shew et al., 2009, 2011), each culture (n = 15) was comprised of a coronal slice of somatosensory rat cortex combined with a slice from the ventral tegmental area, which provides dopaminergic inputs to the cortex for proper network development (Gireesh and Plenz, 2008). The tissue was cultivated directly on the surface of an 8 × 8 grid of electrodes (Fig. 1A, 200 μm interelectrode distance, 30 μm electrode diameter, no corner electrodes, 4 kHz sampling). Recordings were taken between 10 and 20 d in culture allowing several 1 h recordings from each network (47 recordings in total; Stewart and Plenz 2008). The recorded voltages were low-pass filtered at 50 Hz to obtain the LFP, which was shown to correlate with the spiking activity of the local neuronal population near each electrode (Shew et al., 2009).
Spontaneous neuronal burst activity in cortex organotypic cultures grown on integrated planar MEA. A, Light-microscopic image of a neural culture including somatosensory cortex and VTA grown on a 60 channel MEA. Black dots indicate the location of the n = 59 recording electrodes. The ground electrode is outside the imaged region. B, Examples of LFP traces simultaneously recorded during a burst of neural activity (no drugs applied). The spatial arrangement of LFP traces matches that of the electrodes from which they were recorded. Inset: magnified view of three LFP traces. C, Average normalized power spectra for normal (black, no drug, n = 8), disfacilitated (blue, DNQX/AP5, n = 8), and disinhibited conditions (red, PTX, n = 8). The spectra were typically broadband without strong peaks under normal and disfacilitated conditions, while prominent β-oscillations (10–20 Hz) emerged under disinhibited conditions. Error bars indicate ±SEM.
Observed dynamics consisted of bursts of activity like the example shown in Figure 1B, which often spanned many recording sites. We determined the start and end of each burst as well as which sites participated using a threshold to identify times and electrodes with large-amplitude negative LFP deflections (Material and Methods).
We applied pharmacological agents to change the network excitability. Neural synchrony is expected to be sensitive to such changes in excitatory and inhibitory interactions as we demonstrate with our computational model below. Excitatory synaptic transmission was reduced with combined application of the NMDA and AMPA glutamate receptor antagonists AP5 and DNQX. Inhibitory synaptic transmission was reduced with GABAA receptor antagonist PTX. Empirically, we found that the used concentrations of DNQX/AP5 resulted in disfacilitation, i.e., decreased burst size and duration, while the partial disinhibition with PTX increased burst size and duration. As shown in Figure 1C (blue, black), we found that the frequency content of the LFP signals did not show any strong peaks at particular frequencies for the no-drug and DNQX/AP5 conditions. In contrast, β-oscillations became more prominent in disinhibited cultures as reported previously (Fig. 1C, red; Gireesh and Plenz, 2008).
To parameterize these drug effects on spontaneous network dynamics, we used the statistical measure κ, which is based on measured distributions of burst sizes, as developed previously (Shew et al., 2009, 2011). This method takes advantage of the fact that the unperturbed condition typically results in a burst size distribution of the form Pr(s)∼s−1.5, i.e., neuronal avalanches, while the disfacilitated and disinhibited conditions result in exponential and bimodal distributions respectively. Practically, κ ≈ 1 for unperturbed networks (κ = 1.14 ± 0.01; no drug; n = 28), κ < 1 for disfacilitated networks (κ = 0.81 ± 0.01; DNQX/AP5; n = 10) and κ > 1 for disinhibited networks (κ = 1.51 ± 0.01; PTX; n = 9). This relationship is apparent in Figure 2A, where each point represents the average spatial area a of bursts and κ value from a single recording and color represents drug condition. The average burst area was very small and increased gradually for small κ, and began to rise more quickly near κ = 1 and beyond (Fig. 2A). We next quantified the variability of burst area by computing the Shannon entropy H(a) of the burst area distribution. In line with similar previous work (Shew et al., 2011), the entropy peaked near, but slightly higher than κ = 1, demonstrating that the variability of burst area was highest under conditions which favor neuronal avalanches (Fig. 2B). Similar to spatial area, the mean burst duration d began to increase sharply near κ = 1 and the entropy of burst duration H(d) peaked near κ = 1 (Fig. 2C,D).
Burst area and duration have moderate mean and maximum entropy near κ = 1. A, Network dynamics were parameterized using κ. Low κ, κ ≈ 1, and high κ indicate disfacilitated, neuronal avalanche, and disinhibited dynamics respectively. As network activity increased with an increase in κ, the average spatial extent of spontaneous bursts rose slowly for small κ, and more steeply near κ = 1. Each point is the averaged burst area from a 1 h recording. B, The diversity, i.e., Shannon entropy, of burst area was highest near κ = 1. C, The average duration of bursts also rose strongly near κ = 1 and saturated for high κ. D, Entropy of burst duration peaked near κ = 1. Error bars indicate ±SEM.
Neuronal synchrony has moderate mean and maximum entropy near κ = 1
During a burst, the LFP recorded from different electrodes were often, but not always synchronized (Fig. 1B, inset). We hypothesized that experimental conditions which result in neuronal avalanches also result in moderate average synchrony and maximally variable synchrony. To assess synchrony, we first obtained a phase trace for each LFP trace using the Hilbert transform (Fig. 3; Materials and Methods). To visualize the phases of all 59 electrodes versus time, we used dynamic phase histograms in which periods of elevated phase synchrony appear as yellow/red “bundles” of phase trajectories (Fig. 4A, top). This is further clarified in the dynamic phase difference histogram, in which yellow/red pixels near a zero phase difference indicate that many pairs of electrodes have the same phase, i.e., are phase synchronized (Fig. 4A, middle). The temporal changes in phase synchrony in the network were quantified using Kuramoto's order parameter, r (Fig. 3C; Strogatz, 2001; Arenas et al., 2008), which is based on all 59 phases. Periods of high phase synchrony, i.e., low phase difference, corresponded to high values of r (Fig. 4A, bottom). As shown in the examples in Figure 4, the duration of phase synchrony changed with κ from relatively short periods for κ < 1 (Fig. 4B), to moderate and prolonged phase synchrony periods for κ ≈ 1 and κ > 1 (Fig. 4A,C).
Assessing phase synchrony from multisite LFP recordings. A, Example LFP amplitude trace recorded from a single electrode under disinhibited condition. A 100 ms period (red box) is expanded (lower trace). B, Using the Hilbert transform, the LFP amplitude trace is converted to a corresponding phase trace. Top, Expanded phase period. Note that phase near ±π corresponds to negative LFP peaks, while phase near zero coincides with positive LFP peaks (arrows). Lower trace shows corresponding phase for the full time period in A. C, Examples for lack of phase-locking (left), in-phase locking (middle) and combined in-phase and anti-phase locking (right) taken from four sites (colors). Top row, LFP traces. Middle row: phase traces. High phase synchrony appears as times when many sites display a similar phase trace simultaneously (middle), while dispersed phases indicate lack of synchrony (left). Anti-phase locking appears as bundles of in-phase traces separated by π (right; double headed arrow). Note low LFP amplitudes during lack of synchrony. Bottom row, The phases at the time indicated by the black arrow are visualized on a unit circle. The phase synchrony measure, r, is the distance (length of red arrow) between the origin and the center of mass of the phases. Note that anti-phase locking tends to reduce r.
Phase synchrony dynamics for example bursts under different conditions of network excitability. A, Moderate levels of phase synchrony under normal condition. Time of the burst indicated by the black bar. Top, Dynamic phase histogram. Color indicates the number of electrodes with a given phase (vertical axis) at a given time (horizontal axis). Phase bundles marked by yellow/red pixels over time indicate many sites with the same phase, i.e., in-phase locking. Middle, Dynamic phase difference histogram. Color indicates the number of electrode pairs with a given phase difference (vertical axis) at a given time (horizontal axis). Phase differences near 0 indicate in-phase locking. The gray histogram to the right is the time average of the dynamic phase difference histogram. Note dominance of in-phase locking. Bottom, Phase synchrony r. Dashed line indicates r = 0. Phase synchrony r is high during periods of in-phase locking. B, Disfacilitated activity is characterized by brief and weakly synchronized bursts. Same as A, but during an example burst in the presence of DNQX/AP5. C, Disinhibited activity is characterized by strongly synchronized long-lasting bursts, anti-phase locking (phase differences near π), and an underlying prominent ∼20 Hz β-oscillation. Same as A, but in the presence of PTX. Burst duration extends beyond period shown. Quantification of anti-phase locking, APL, is illustrated in time-averaged dynamic phase difference histogram (right).
To quantify these changes in synchrony with κ, we first calculated the average network synchrony per burst, S̄N. This was done by integrating r for the duration of each burst and averaging over all bursts for each experiment. We found that, as the network was changed from disfacilitated to disinhibited, S̄N increased slowly at first, began to rise more sharply near κ = 1, and reached a plateau for large κ (Fig. 5A). Next we quantified the diversity of network synchrony, H(SN). We found that H(SN) is maximized near, but just above κ = 1 (Fig. 5B) indicating that networks with neuronal avalanches achieve the most diverse repertoire of network synchronization during bursts.
Neural synchrony attains moderate mean and maximum entropy near κ = 1. A, Mean network synchrony showing a slow rise for low κ, a sharp increase near, just above κ = 1, and saturation for high κ. B, The entropy of network synchrony showing a peak close to κ = 1. C, Mean instantaneous network synchrony showing a gradual rise with κ indicating that the increase in burst duration does not fully explain the rising trend in A. Broken line is a control computed during duration-matched non-burst periods. D, Entropy of instantaneous network synchrony is maximal near κ = 1. E, Mean instantaneous burst synchrony does not change with κ. Thus, burst area and duration are sufficient to explain rising trend in mean synchrony. Broken line: control, computed for electrodes which do not participate in bursts. F, The peak in entropy persists for instantaneous burst synchrony.
The rise in network synchrony as a function of κ, a priori, could depend on both the burst duration and spatial burst area. Therefore, we next explored the degree to which the sigmoidal rise in network synchrony could be explained by the similar κ-dependent changes in the spatiotemporal boundaries of bursts presented in Figure 2. First, for each burst, we normalized network synchrony by burst duration to obtain the average instantaneous network synchrony, S̄IN (Fig. 5C,D) and the corresponding entropy, H(SIN). We found that both the rising trend in synchrony and the peak in entropy at κ ≈ 1 persisted, indicating that burst duration alone is insufficient to explain the trends shown in Figure 5, A and B. Next, to determine the influence of burst area, we recomputed the instantaneous network synchrony for each burst, but included only those electrodes that participated in the burst, which we call instantaneous burst synchrony, SIB (Fig. 5E,F). We note that SIB was also corrected to account for the different expected noise levels for different numbers of electrodes, so that bursts covering different areas were fairly compared. For comparison, the instantaneous out-of-burst synchrony among those electrodes that did not participate in the burst was significantly lower (Fig. 5E, broken line). Because S̄IB did not show a significant rising trend versus κ, we concluded that the trend in Figure 5A is primarily determined by a combination of κ-dependent changes in burst area and duration. This result indicates that mean synchrony within the spatiotemporal boundaries of bursts does not significantly change with network excitability. However, the entropy of instantaneous burst synchrony, H(SIB) continued to show a peak. This indicates that entropy of synchrony arises not only from variability of the spatiotemporal boundaries of the bursts, but also from intrinsic, within-burst variability that is maximum near κ = 1.
Importantly, our main findings in Figure 5 were insensitive to changes in (1) the spatial extent of the MEA, (2) the spatial resolution of the MEA, and (3) the frequency band of LFP. The first two of these controls was performed by repeating the analysis using subsets of electrodes from the original recordings. To test for robustness to changes in spatial extent, we used a 4 × 4 set of electrodes near the center of the MEA (Fig. 6, red). The spatial area of this subset is reduced by a factor of 4. To test for robustness to changes in spatial resolution, i.e., interelectrode distance, we used a 4 × 4 subset of electrodes with 400 μm interelectrode distance (Fig. 6, blue). Thus, the spatial resolution is halved compared with the original recording. In both cases, the trends presented in Figure 5 were largely unchanged. This suggests that our study should be repeatable with other MEA systems and our findings reflect cortical dynamics which span a range of spatial scales.
Synchrony as function of κ is independent of spatial extent, resolution of MEA, and frequency band of LFP. The data shown and described in Figure 5 were recomputed using two different arrangements of MEA electrodes and three different frequency bands. A, A 4 × 4 compact group of electrodes covering a 800 μm × 800 μm recording area near the center of the array (red). B, A 4 × 4 sparse group of electrodes covering the full 1600 μm × 1600 μm area, but with half the spatial resolution (blue). C, Frequency bands included 4–12 Hz (green), 12–25 Hz (gold), and 25–50 Hz (orange). D–I, Same figure legend as in Figure 5A–F. Note that main conclusions were unchanged, although average network synchrony tended to decrease for higher frequency bands.
Robustness to changes in frequency band was demonstrated by repeating our analysis for three bands: 4–12, 12–25, and 25–50 Hz. As shown in Figure 6, qualitative changes in phase synchrony versus κ were unchanged for these three bands compared with our original broadband analysis. However, we note that the quantitative levels of synchrony were lower for higher frequency bands in the high κ regime.
Anti-phase locking in disinhibited cortical networks with κ > 1
At high κ, anti-phase locking (Fig. 3C) was found for some bursts, which was visible as a tendency for π phase differences (Fig. 4C). We found that anti-phase locking was ∼10 times more likely in high κ experiments compared with those with low or near κ = 1 (ANOVA, p < 0.01). For the 8 experiments with highest κ, we found APL = 6.3 ± 1.8 (mean ± SEM), where APL quantifies the degree of anti-phase locking, whereas APL = 0.9 ± 0.3 and APL = 0.4 ± 0.1 for the 8 experiments near κ = 1 and low κ respectively (Fig. 4C; Materials and Methods). The fraction fAPL of anti-phase-locked sites, which extended from deep to superficial layers, was fAPL = 0.64 ± 0.02 for the high κ experiments, which is 25% higher than expected by chance and significantly higher than experiments with low κ (fAPL = 0.51 ± 0.001) and κ near 1 respectively (fAPL = 0.52 ± 0.01; ANOVA, p < 0.01). Although such anti-phase locking in general will reduce the magnitude of r, the high κ conditions still resulted in the highest average r due to the fact that anti-phase-locked groups were always small compared with the phase-locked groups.
Maximum entropy and onset of synchrony near κ = 1 in computational model
To gain further insight into our experimental findings, we next performed a similar investigation in a network-level computational model. Network-level phase synchrony has been investigated extensively in previous theoretical and model studies using coupled oscillators (Strogatz, 2001; Arenas et al., 2008). Naively, one might interpret the signal recorded from one electrode as one theoretical oscillator. However, in reality, each electrode records the collective activity of a population of E and I neurons. Thus, our experimental pharmacological E–I manipulations effect not only the coupling between electrodes, but also the coupling among the local neurons measured by a single electrode. Since the intrinsic properties of a theoretical oscillator do not depend on coupling, the interpretation of one electrode as one oscillator is problematic. Therefore, we developed a computational model with the aim of obtaining results which were more directly comparable to our experiments, while keeping the model as simple as possible.
The computational model was comprised of an all-to-all connected network of binary neurons, 80% excitatory and 20% inhibitory. The model neurons were divided into 59 groups meant to represent the populations of neurons measured by our 59 electrodes in the experiments (Fig. 7A). Different network sizes were studied with between 30 and 80 neurons per group. Probabilistic integrate-and-fire dynamics were implemented (Materials and Methods). We modeled bursts like those observed in our experiments one at a time, each starting with a single initially active excitatory neuron and continuing until activity ceased or a time limit was reached. The spike activity of each group was summed to model an LFP-like population signal (Fig. 7B,C). This is justified based on previous experiments showing that the size of a burst based on LFP measurements is proportional to the number of spikes recorded from the underlying network during the burst (Shew et al., 2009). These group signals were then subjected to precisely the same analysis as performed with the experimental data. Figure 7, D and E, shows example phase traces, a dynamic phase histogram, and a synchrony r trace. We quantified how burst area, duration, and phase synchrony depend on E–I conditions, as revealed by κ.
Phase synchrony and neuronal avalanches in a network-level computational model of E–I neurons. A, A recurrently connected population of 80% excitatory and 20% inhibitory neurons were divided into 59 groups with Q neurons in each group (Q = 30, 40, 50, 60, 70, 80 were tested). Each group models the neurons nearby one recording site in the experiment. B, Spike rasters from two groups. (blue, E; red, I). C, Population signals corresponding to the two group rasters shown in B. The spikes generated by a group at each time step were summed (light colored traces) and bandpass filtered (dark colored traces) to generate population signals. D, Phase of these filtered population signals was obtained using the Hilbert transform. E, An example dynamic phase histogram and synchrony trace. Phase synchrony among the 59 groups was analyzed exactly as in the experiments. F, Bursts of activity were modeled one at a time. Burst size distributions were used to compute κ. Neuronal avalanches (κ ≈ 1) were observed only when E and I were balanced such that, on average, the spike count of the population did not grow or decay over consecutive time steps. Bimodal size distributions (κ > 1) occurred when inhibition was suppressed. And a lack of large-sized bursts marked the disfacilitated distribution (κ < 1).
Disinhibited or disfacilitated conditions were modeled by weakening all efferent connections from inhibitory or excitatory neurons, respectively. The baseline, normal condition of the computational model was chosen such that activity in the network propagated in a balanced manner. More specifically, we established the magnitude of connections between neurons such that a single active excitatory neuron was expected, on average, to excite exactly one more neuron in the following time step. Disinhibition and disfacilitation disrupted this balance. Disinhibition made it likely that more than one neuron would be excited and resulted in growth (followed eventually by saturation) of the number of active neurons over time, on average. Disfacilitation resulted in decay of the number of active neurons over time.
The distribution of burst sizes under the balanced normal condition was a power law with exponent near −1.5 (Fig. 7F), i.e., the computational model generated neuronal avalanches and κ ≈ 1. The burst size was defined as the number of spikes that occurred during the burst. When disinhibited, our computational model produced bursts with a bimodal probability distribution, i.e., κ > 1, and when disfacilitated the distribution was closer to an exponential, κ < 1, in good agreement with our experiments (Fig. 7F).
We found that as the computational model was tuned from low to high κ, the characteristics of bursts matched those observed in our experiments. Most importantly, average burst area, duration and network synchrony all displayed a sharp rise near the critical balanced regime identified by κ = 1 (Fig. 8A–C, left). The entropy of these quantities reached a peak near, but just above κ = 1, as observed in the experiments (Fig. 8A–C, right). Instantaneous network synchrony (Fig. 8D) and instantaneous burst synchrony (Fig. 8E) versus κ also matched the experiments well. However, there were quantitative differences between the computational model and the experiments. First, the onset of average quantities and peak entropy for the model was typically shifted to slightly larger κ values than those observed experimentally. Second, the range of obtainable κ values was lower in the model—even with inhibition completely suppressed, the model did not reach κ = 1.6 as found for the most extreme experiments. Finally, we did not observe anti-phase locking in our model. Such differences could arise from many possible simplifications made in the computational model, including all-to-all connectivity, no large-scale network structure like cortical layers, and binary neural state. However, we emphasize that our main conclusions hold for both the experiments and the computational model; namely, peak entropy and moderate average synchrony occurred under E–I conditions which result in neuronal avalanches.
Onset of synchrony and peak entropy coincide near the critical balance of E–I in computational model. All of the experimental results presented above were well matched by the computational model. A, Average burst area (left) increased around κ ≈ 1 where entropy was highest (right). B, Average burst duration (left) increased near κ = 1 where entropy was highest (right). C–E, Same as captions from Figure 5A–F. The largest difference between the experiments and the model was the shift in peak entropy to higher κ values for SIN in the model. The dashed lines delineate the variation (minimum to maximum) observed over 60 different simulated networks (10 networks at each of 6 sizes).
Discussion
The ability of a network of cortical neurons to synchronize its dynamics depends on recurrent interactions among excitatory and inhibitory components of the network. Here we studied spontaneously emerging network-level synchrony over a range of E–I conditions in rat cortex slice cultures and a computational model. We demonstrated the existence of a critical level of network excitability identified by κ ≈ 1. Below criticality, synchrony is minimal; just above criticality, significant synchrony emerges; well above criticality, synchrony reaches intense levels. Close to criticality, variability (entropy) of synchrony was highest and neuronal avalanches were observed.
Below, we first discuss potential implications of variable synchrony. We then relate our high κ experiments to epilepsy models. Next, we discuss our findings of anti-phase locking. Finally, we highlight similarities between our results and theoretical predictions from critical phenomena, which suggest that the κ ≈ 1 condition corresponds to the critical point of a phase transition.
Variability of synchrony
One implication of our results derives from the fact that either abnormally low or high spontaneous synchrony is associated with cortical dysfunction. Our observations of a sigmoidal rise in average synchrony as a function of network excitability suggests that a limited range of E–I conditions (≈1 < κ < 1.25) will confer the cortex with a healthy, moderate range of mean synchrony. This same range of excitability resulted in the highest variability of synchrony, suggesting that extremely variable synchrony may be unavoidable if the cortex must operate with moderate mean synchrony.
The implications of highly variable spontaneous synchrony for cortical information processing are not well understood. One traditional view is that such “correlated noise” impedes reliable coding of input from external stimuli. If this view is correct, our results highlight a competition: moderate average synchrony comes at the cost of high variability. Indeed, fluctuations of ongoing activity influence the processing of input (Arieli et al., 1996; Poulet and Petersen, 2008; Kelly et al., 2010). However, a second view suggests this influence is not necessarily detrimental, particularly when these fluctuations occur with balanced E–I (Cafaro and Rieke, 2010). A third possible purpose of spontaneous activity is that it supports an “internal model” of previous experiences (Berkes et al., 2011), which is consistent with many observations of similarities between spontaneous activity and previous stimulus-evoked activity (Tsodyks et al., 1999; Ji and Wilson, 2007; Han et al., 2008; Luczak et al., 2009). In this view, maintaining cortical E–I conditions with high entropy of synchrony, i.e., avoiding high or low κ, may be crucial. Indeed, an internal model which accurately reflects diverse and variable previous experiences may require a diverse and variable (high entropy) repertoire of spontaneous synchrony.
We note that the peak in entropy of synchrony we observed experimentally persisted even after the variability in burst duration and spatial area were accounted for (Figs. 5F, 6I). This was also confirmed in our computational model (Fig. 8D,E, right). This implies that even if one hand-selected a subset of bursts which all have the same duration and area, their network synchrony would vary from burst to burst. This variability would peak near κ = 1. This suggests that the peak in entropy results in part from variability at the single-electrode level. This is understandable considering that, for both our experiment and our computational model, each electrode samples a subnetwork consisting of tens to hundreds of E and I neurons. Our manipulations of E–I conditions must effect synchrony within such single-site micro-circuits. Thus, the same E–I conditions which result in variable synchrony at the large scale (between electrodes) are also likely to result in variability synchrony at the micro-circuit scale, which may explain the peak in entropy of instantaneous burst synchrony.
κ > 1 identifies a pathological state of synchrony typical for epilepsy
Cortical disinhibition is a model for epilepsy (Luhmann et al., 1995; Prince et al., 2009). GABAA antagonists readily induce epileptiform seizures in isolated cortex preparations (Gutnick et al., 1982) that initiate in and involve deep layer pyramidal neurons (Connors, 1984; Pinto et al., 2005). Provided layers 5/6 are intact (Telfeian and Connors, 1998), ∼10–20% reduction in GABAA receptor function suffices to induce seizures (Chagnac-Amitai and Connors, 1989), in line with the relatively low dose of 5 μm PTX concentration used in the current study, which is far from the 50 μm required to block the GABAA receptor. Our observation of spatially extended bursts including deep layers for large κ is in accordance with these studies. Moreover, for high κ, we observed increased β-oscillations (Fig. 1; Gireesh and Plenz, 2008), which are known to involve layer 5/6 pyramidal neurons (Yamawaki et al., 2008). Importantly, the involvement and emergence of β-oscillations in epilepsy has been well documented. Ten to 15 Hz oscillations emerge spontaneously in isolated cortex preparations using the epilepsy model of low extracellular magnesium (Flint and Connors, 1996) and epileptic seizures recorded in animals during slow-wave sleep have elevated β-power (Timofeev and Steriade, 2004). In epileptic patients, long-range temporal correlations are enhanced in the β-band of intracranial EEG recorded during seizure-free epochs (Monto et al., 2007), and entrainment of photosensitive seizures is particularly effective with flashes between 15 and 20 Hz (Parra et al., 2005). We conclude that the regime of κ > 1 identifies a pathological state of relatively large-scale and robust phase synchrony across cortical layers with dominant oscillation frequencies similar to those found in epilepsy.
Anti-phase synchrony for κ > 1
One aspect of our experiments which is unexplained by our computational model is the anti-phase locking at high κ. Anti-phase LFP oscillations could arise from several possible mechanisms. First, in-phase synchronized spiking in a population of cortical pyramidal neurons with aligned long apical dendrites can produce such LFP signals (Chrobak and Buzsáki, 1998). Accordingly, the incidence of anti-phase locking is expected to increase with the involvement of deep layer neurons in generating β-oscillations at high κ. Second, 2-oscillator models reveal anti-phase locking to arise from excitatory coupling (Kopell and Somers, 1995; Ermentrout and Kleinfeld, 2001; Neltner and Hansel, 2001) such as mutual layer-V pyramidal neuron interactions. Third, in line with our findings at high κ, i.e., disinhibition, network simulations demonstrate that inhibition can suppress anti-phase synchrony (Kanamaru, 2006). The concurrence of slow β-oscillations and anti-phase locking in our experiments is in line with anti-phase locking in human motor coordination at lower frequencies as predicted by the HKB-model (Haken et al., 1985; Schöner et al., 1986).
κ ≈ 1 suggests a critical phase transition
Neuronal avalanches were originally named (Beggs and Plenz, 2003) in reference to critical phenomena—a topic which is traditionally concerned with the statistical physics of phase transitions (Hinrichsen, 2006), but has also been applied to understand network-level neural dynamics (Levina et al., 2007; Buice and Cowan, 2009; Millman et al., 2010). The most prominent finding of our experiments, which was predicted by theory of critical phenomena, is the power-law probability distribution of burst sizes, which we observed here when no drugs were applied in our experiments and when E–I was balanced in the computational model. We found a power-law exponent near −1.5, which suggests that cortical networks belong to the class of non-equilibrium critical phenomena called “directed percolation” (Hinrichsen, 2006). Originally devised for the study of a fluid percolating through a porous material, this theoretical framework may also be interpreted as neural activation propagating through a network. Models of this type have successfully explained measurements of cortical dynamic range (Kinouchi and Copelli, 2006; Shew et al., 2009; Larremore et al., 2011) and information transfer (Beggs and Plenz, 2003; Shew et al., 2011). Our own computational model is very closely related to these, but has been modified to include both E and I neurons. Previous work included only E neurons. Inhibitory neurons make our computational model more similar to real cortical networks and facilitate comparison with experiments. However, the implications of inhibitory neurons for the theory of critical phenomena are not yet well understood. Can such an E–I network be truly critical, or perhaps only “near” critical? The observed onset of synchrony together with peak entropy and the power law burst size distribution strongly suggest that our computational model (and by association, our experiments) undergoes a critical phase transition near κ = 1. However, careful theoretical work which accounts for inhibitory neurons will be required for more direct support of this hypothesis.
In summary, our experiments are in agreement with predictions from our computational model, suggesting that, when E–I is unperturbed by drugs, the cortex tends to operate with an intermediate excitability, balanced such that population activity propagates without growth or decay on average. This balanced excitability marks the boundary between a weak synchrony regime and an intense synchrony regime. Neuronal avalanches and maximal variability of synchrony also lie at the boundary between these regimes. Normal, healthy cortex seems to operate with moderate synchrony on average. Our results suggest that to achieve such moderate synchrony, the cortical E–I must be regulated to the balanced excitability indicated by neuronal avalanches (κ ≈ 1).
Footnotes
This work was supported by the Intramural Research Program of the National Institute of Mental Health. R.R. acknowledges support from Department of Defense Multidisciplinary University Research Initiative Grant ONR N000140710734 to the University of Maryland. We thank Craig Stewart for help preparing the cultures and Dr. Shan Yu for helpful discussions.
- Correspondence should be addressed to Dr. Dietmar Plenz, Section on Critical Brain Dynamics, Laboratory of Systems Neuroscience, National Institute of Mental Health, Bethesda, MD 20892. plenzd{at}mail.nih.gov
