Abstract
Spontaneous neuronal activity is a ubiquitous feature of cortex. Its spatiotemporal organization reflects past input and modulates future network output. Here we study whether a particular type of spontaneous activity is generated by a network that is optimized for input processing. Neuronal avalanches are a type of spontaneous activity observed in superficial cortical layers in vitro and in vivo with statistical properties expected from a network operating at “criticality.” Theory predicts that criticality and, therefore, neuronal avalanches are optimal for input processing, but until now, this has not been tested in experiments. Here, we use cortex slice cultures grown on planar microelectrode arrays to demonstrate that cortical networks that generate neuronal avalanches benefit from a maximized dynamic range, i.e., the ability to respond to the greatest range of stimuli. By changing the ratio of excitation and inhibition in the cultures, we derive a network tuning curve for stimulus processing as a function of distance from criticality in agreement with predictions from our simulations. Our findings suggest that in the cortex, (1) balanced excitation and inhibition establishes criticality, which maximizes the range of inputs that can be processed, and (2) spontaneous activity and input processing are unified in the context of critical phenomena.
Introduction
The cortex is spontaneously active even in the absence of any obvious stimulus or motor output. Increasingly, evidence shows that such ongoing activity is intricately linked to stimulus-evoked activity. For example, orientation maps constructed from ongoing neuronal activity in the anesthetized cat match those based on visual responses (Tsodyks et al., 1999; Kenet et al., 2003). Spatiotemporal correlations of spikes in the visual cortex are similar when the awake animal is simply sitting in darkness or observing natural scenes (Fiser et al., 2004). Likewise, population responses to auditory and somatosensory stimuli fall within the repertoire of observed spontaneous events (Luczak et al., 2009). Moment to moment, ongoing activity contributes to the large variability observed in stimulus responses (Arieli et al., 1996; Azouz and Gray, 1999; Kisley and Gerstein, 1999), while being only weakly modulated by stimulus presentation (Fiser et al., 2004). On longer timescales, the organization of spontaneous activity is thought to reflect past inputs and influence future network responses (Ohl et al., 2001; Yao et al., 2007). Such interplay between spontaneous and stimulus-evoked activity raises the question whether there is a particular type of ongoing activity that maintains optimized stimulus processing in the network.
Here we focus on neuronal avalanches, a type of spontaneous activity observed in superficial layers of cortex in vivo and in vitro (Beggs and Plenz, 2003; Plenz and Thiagarajan, 2007; Gireesh and Plenz, 2008; Petermann et al., 2009). Neuronal avalanches consist of bursts of elevated population activity, correlated in space and time, that are distinguished by a particular statistical character: activity clusters of size s occur with probability P(s) ∼ sα, i.e., a power law with exponent α = −1.5. Neuronal avalanches have several additional key properties: (1) they arise during development when superficial layers form in vitro and in vivo (Gireesh and Plenz, 2008), (2) they are homeostatically maintained for weeks in isolated cortex without any input (Stewart and Plenz, 2008), (3) they constitute the dominant form of ongoing cortical activity in the awake monkey (Petermann et al., 2009), and (4) their pharmacological regulation is characterized by an inverted-U profile of NMDA/dopamine D1 receptor interaction and intact fast inhibitory transmission (Beggs and Plenz, 2003; Stewart and Plenz, 2006; Gireesh and Plenz, 2008).
Neuronal avalanches are similar to the dynamics of other systems poised at the boundary of order and disorder; more precisely, we refer to systems operating at criticality (Stanley, 1971; Bak and Paczuski, 1995; Jensen, 1998). Importantly, simulations predict that at criticality, neuronal networks optimize several aspects of information processing, including (1) the range of stimulus intensities that can be processed, i.e., dynamic range (Kinouchi and Copelli, 2006), and (2) the amount of information that can be transferred (Beggs and Plenz, 2003; Tanaka et al., 2009). Until now, experimental support of these predictions has been lacking. Here, we demonstrate that in vitro cortical networks have maximum dynamic range when spontaneous activity takes the form of neuronal avalanches. By systematically changing excitation and inhibition, we obtain a tuning curve for stimulus processing in cortical networks, with peak performance found under balanced conditions that generate neuronal avalanche activity.
Materials and Methods
Organotypic cultures on microelectrode arrays.
Coronal slices from rat somatosensory cortex (350 μm thick, postnatal day 0–2; Sprague Dawley) and the midbrain (VTA; 500 μm thick) were cut and cultured on a poly-d-lysine-coated 8 × 8 microelectrode array (MEA) (Multi Channel Systems; 30 μm electrode diameter; 200 μm interelectrode distance). In this organotypic slice coculture, the development of deep and superficial cortical layers (Götz and Bolz, 1992; Plenz and Kitai, 1996) as well as neuronal avalanche activity parallels that observed in vivo (Gireesh and Plenz, 2008). In short, a sterile chamber attached to the MEA allowed for repeated recording from cultures for weeks. After plasma/thrombin-based adhesion of the tissue to the MEA, standard culture medium was added (600 μl, 50% basal medium, 25% HBSS, 25% horse serum; Sigma-Aldrich). MEAs were then affixed to a slowly rocking tray inside a custom-built incubator (±65° angle, 0.005 Hz frequency, 35.5 ± 0.5°C) (Stewart and Plenz, 2008).
Spontaneous activity.
Using a recording head stage inside the incubator (MEA1060 w/blanking circuit; ×1200 gain; bandwidth 1–3000 Hz; 12 bit A/D; range 0–4096 mV; Multi Channel Systems), local field potential (LFP; 4 kHz sampling rate; reference electrode in bath) was obtained from 1 h of continuous recordings of extracellular activity (low-pass, 100 Hz, phase-neutral). To establish a correlation between LFP and neuronal spiking activity, in n = 5 cultures, extracellular activity was recorded for 15 min at 25 kHz. In addition to extracting LFP, the extracellular signal was filtered in the frequency band 300–3000 Hz and <78 single units were identified per culture using threshold detection and PCA-based spike sorting (Offline Sorter; Plexon).
Stimulus-evoked activity.
Immediately following each 1 h recording of spontaneous activity, stimulus-evoked activity was measured. Stimuli were applied at 5 s intervals at one electrode located approximately at the center of the culture, in superficial cortical layers. Stimuli were current-controlled, single shocks with bipolar square waveform: 50 μs with amplitude −S followed by 100 μs with amplitude +S/2, where 6 < S < 200 μA. We tested one set of stimulus amplitudes with fine resolution (S = 10–200 μA in steps of 10 μA), and another with coarser resolution (S = 6, 12, 24, 50, 65, 80, 100, 150, 200 μA). Results were similar for the two protocols. Each stimulus level was repeated 40 times in pseudorandomized order resulting in a total recording duration of 2000 (coarse) or 4000 (fine) s. Each stimulus-evoked response was recorded using all electrodes except for the stimulation electrode during 500 ms following stimulation. A blanking circuit disconnected the recording amplifiers during stimulation, significantly reducing stimulus artifacts (Multi Channel Systems). Sample rate and filtering was identical to that used for spontaneous activity recordings.
Pharmacology.
Bath application of antagonists of fast glutamatergic or GABAergic synaptic transmission was used to change ratios of excitation to inhibition (E/I). The normal (no-drug) followed by a drug condition was studied within 3 h to minimize nonstationarities during development. Stock solutions were prepared for the GABAA receptor antagonist picrotoxin (PTX), the NMDA receptor antagonist (2R)-amino-5-phosphonovaleric acid (AP5), and the AMPA receptor antagonist 6,7-dinitroquinoxaline-2,3-dione (DNQX). Six microliters of these stock solutions were added to 600 μl of culture medium to reach the following working concentrations (in μm): 5 PTX, 20 AP5, 10 AP5 + 0.5 DNQX, and 20 AP5 + 1 DNQX. After recording, the drug medium was replaced with 300 μl of drug-free conditioned medium (collected from the same culture the previous day) mixed with 300 μl of fresh, unconditioned medium. Most cultures recovered to criticality within ∼24 h.
Spontaneous cluster size and response to stimulus.
For each electrode, we identified negative peaks in the LFP (nLFPs) that were more negative than −4 SDs of the electrode noise. We then identified a cluster of nLFPs on the array as a group of consecutive nLFPs each separated by less than a time τ (Beggs and Plenz, 2003). The threshold τ was chosen to be greater than the short timescale of interpeak intervals within a cluster, but less than the longer timescale of intercluster quiescent periods (τ = 86 ± 71 ms for all cultures; see also supplemental material, available at www.jneurosci.org). Results were robust for a large range in the choice of τ (data not shown). The size s of a cluster was quantified as the absolute sum of all nLFP amplitudes within a cluster. Similarly, the size R of an evoked response was quantified as the absolute sum of nLFPs within 500 ms following a stimulus.
Definition of κ.
For neuronal avalanches, the probability density function (PDF) of cluster size s follows a power law with slope α = −3/2 (Beggs and Plenz, 2003) (see Fig. 2A). Thus, the corresponding cumulative density function (CDF) for cluster sizes, FNA(β), which specifies the fraction of measured cluster sizes s < β, is a −1/2 power-law function, FNA(β) = (1 −
where βk are m = 10 burst sizes logarithmically spaced between the minimum and maximum observed burst size. Using CDFs rather than PDFs to calculate κ avoids the sensitivity to binning in constructing a PDF. Compared to other nonparametric comparisons of CDFs, e.g., Kolmogorov–Smirnov and Kuiper's test, as well as other methods, κ more accurately measures deviation from neuronal avalanches (see supplemental material, available at www.jneurosci.org).
Dynamic range.
After measuring responses to a range of stimulus amplitudes, we used the response curve, R(S), to compute dynamic range,
where Smax and Smin are the stimulation values leading to 90% and 10% of the range of R, respectively.
Model.
The model consisted of N all-to-all coupled, binary-state neurons (N = 250, 500, 1000) and the following dynamical rules: If neuron j spiked at time t (i.e., sj(t) = 1), then postsynaptic neuron i will spike at time t + 1 with probability pij. As such, the pij are N2 numbers representing the synaptic coupling strengths between each pair of neurons. The pij are asymmetric pij ≠ pji, positive, time-independent, uniformly distributed random numbers with mean and SD of order N−1. If a set of neurons J(t) spikes at time t, then the probability that neuron i fires at time t + 1 is exactly piJ(t) = 1 − ∏j∈J(t)(1 − pij). To implement the probabilistic nature and variability of unitary synaptic efficacy, neuron i actually fires at time t + 1 only if piJ(t) > ζ(t), where ζ(t) is a random number from a uniform distribution on [0,1],
where θ[x] is the unit step function. Like our experiments, we explore a range of network excitability by tuning the mean value of pij from 0.75/N to 1.25/N in steps of 0.05/N by scaling all pij by a constant. For such small mean pij, the model reduces to probabilistic integrate-and-fire, i.e., piJ ≈ Σj∈J(t)pij to order N−2 accuracy. If the mean pij is exactly N−1, then n spikes at time t will, on average, excite n postsynaptic spikes at time t + 1, which constitutes criticality in our model (Beggs and Plenz, 2003; Kinouchi and Copelli, 2006). When mean pij is larger than or less than N−1, the system is supercritical or subcritical, respectively. We define the control parameter of the model σ 
N−1ΣiΣipij. In the context of dynamics, σ reflects the average ratio of spiking descendants to spiking ancestors in consecutive time steps. At criticality, σ = 1; the coupling strengths are balanced such that, on average, the number of active sites neither grows nor decays with time (note that the instantaneous activity level fluctuates greatly). To obtain response as a function of stimulus in the model, we simulated increasing stimulus amplitude S by increasing the number of initially activated neurons (S = 1, 2, 4, 16, 32, 64, 128 initially active neurons). Finally, we note that our model is very similar to (N − 1)-dimensional directed percolation (Buice and Cowan, 2007). Therefore, at high dimension (N > 5) and weak coupling, it is expected that the model behaves as a branching process, where σ is the branching parameter and the −3/2 power law is predicted at criticality.
To test for statistical differences between groups, a one-way ANOVA followed by a Tukey post hoc test was used.
Results
Cortex–VTA cocultures from rat (n = 16), which closely parallel in vivo differentiation and maturation of cortical superficial layers (Gireesh and Plenz, 2008), were grown on 8 × 8 integrated planar microelectrode arrays (Fig. 1A). LFPs were recorded after superficial layer differentiation (>10 d in vitro, DIV) and analyzed to extract spatiotemporal clusters of nLFPs (n = 47 experiments). Extracellular unit activity recorded simultaneously with the LFP revealed that sizes of nLFP clusters correlated with the level of suprathreshold neuronal activity in the network (Fig. 1B) (R = 0.84 ± 0.13, mean ± SD; n = 5 cultures). For each experimental condition, we first measured spontaneous activity (Fig. 1C) and quantified the deviation of the observed spontaneous network dynamics from neuronal avalanche dynamics by calculating κ (Fig. 2A) (see Materials and Methods). In a second step, we measured the input/output dynamic range Δ of the cultured network based on its response to a range of stimulus amplitudes (Figs. 1D, 3A). These measurements were performed under normal conditions and repeated after changing the ratio of excitation and inhibition through bath application of the antagonists PTX or AP5/DNQX.
Measuring spontaneous and stimulus-evoked activity from cortical networks. A, Light-microscopic image of a somatosensory cortex and dopaminergic midbrain region (VTA) coronal slice cultured on a 60-channel microelectrode array. Yellow dot, Stimulation site. Black dots, Recording sites. B, Number of extracellular spikes correlates with the size of simultaneously recorded nLFP burst (R = 0.84 ± 0.13; n = 1). Each point represents total number of spikes versus the corresponding spontaneous nLFP burst size. C, Example recordings of spontaneous LFP fluctuations (left) and nLFP rasters (right) for three drug conditions (top, AP5/DNQX; middle, no drug; bottom, PTX). D, Examples of LFP evoked by 70 μA stimulus (left) and rasters recorded during the application of four stimuli of amplitudes 50, 40, 90, 150 μA (yellow line, stimulus time) (right) for three drug conditions. For both spontaneous (C) and stimulus-evoked (D) activity AP5/DNQX (PTX) typically results in reduced (increased) amplitude LFP events with lesser (greater) spatial extent. In C and D, black dots on the LFP traces indicate nLFP events, raster point color indicates nLFP amplitude, and all calibration bars (left) represent 50 μV, 100 ms.
Change in the ratio of excitation/inhibition moves cortical networks away from criticality. A, Left, PDFs of spontaneous cluster sizes for normal (no-drug, black), disinhibited (PTX, red), and hypoexcitable (AP5/DNQX, blue) cultures. Broken line, −3/2 power law. Cluster size s is the sum of nLFP peak amplitudes within the cluster; P(s) is the probability of observing a cluster of size s. Right, Corresponding CDFs and quantification of the network state using κ, which measures deviation from a −1/2 power law CDF (broken line). Vertical gray lines, The 10 distances summed to compute κ, shown for one example PTX condition (red). B, Simulated cluster size PDFs (left) and corresponding CDFs (right) for different values of the model control parameter σ. C, Summary statistics of average κ values for normal, hypoexcitable, and disinhibited conditions (*p < 0.05 from normal). D, In simulations, κ accurately estimates σ. Broken line, κ = σ. Colored dots, Examples shown in B.
Stimulus–response curves and dynamic range Δ. A, Experimental response R evoked by current stimulation of amplitude S for three example cultures with different κ values. Orange arrows, Range from Smin to Smax; length is proportional to Δ. Note that Δ is largest for κ ≅ 1. B, Model response evoked by different numbers of initially activated sites; Δ is largest for σ ≅ 1. Like the experiment, each point is calculated from 40 stimuli. Error bars, 1 SE. C, Experimental summary statistics for κ under different pharmacological conditions (*p < 0.05 from normal). D, Simulation summary statistics for κ comparing different ranges of σ (*p < 0.05 from σ ≅ 1).
Quantifiying the cortical network state based on κ
Figure 2A (left) shows experimental cluster size PDFs obtained from three cultures under normal, unperturbed conditions and in the presence of PTX or AP5/DNQX, respectively. Under normal conditions, cultures revealed a PDF close to −3/2 power law as predicted for neuronal avalanches (Fig. 2A, black). In the presence of PTX, however, the PDF is bimodal, revealing a high likelihood for small and large activity clusters, but a decreased probability of medium-sized clusters (Fig. 2A, red). In contrast, bath application of AP5/DNQX reduces large clusters resulting in mostly small clusters (Fig. 2A, blue). These differences in PDFs are robustly assessed using the corresponding CDFs (Fig. 2A, right). Reducing excitation results in a steep early rise of the CDF, while reducing inhibition results in a delayed rise of the CDF. κ robustly quantifies these observations using the difference between a measured CDF of cluster sizes and the theoretically expected reference CDF for neuronal avalanches (Fig. 2A, right, gray lines). As summarized in Figure 2C, κ ≅ 1 under normal conditions (κ = 1.14 ± 0.01, mean ± SE; n = 28), κ < 1 when excitation is reduced (κ = 0.81 ± 0.01; n = 10) and κ > 1 when inhibition is reduced (κ = 1.51 ± 0.01; n = 9; F(2,44) = 82.7; p < 0.05 for PTX and AP5/DNQX from normal).
This experimental strategy was paralleled using a network-level computational model of binary, integrate-and-fire neurons, in which changes in the excitation/inhibition ratio (E/I) were mimicked by tuning the parameter σ (see Materials and Methods). For σ < 1, a neuron triggers activity in less than one neuron, on average, resulting in a hypoexcitable state. Conversely, for σ > 1, one neuron excites on average more than one neuron in the near future, resulting in a hyperexcitable condition. Accordingly, for σ = 1, propagation of activity is balanced, as was found experimentally for neuronal avalanches (Beggs and Plenz, 2003; Stewart and Plenz, 2008). We simulated “spontaneous” activity clusters by activating a single randomly chosen neuron and monitoring the ensuing until activity ceased or 500 time steps were executed. The total number of spikes in a cluster was taken as the cluster size. One thousand clusters were simulated at each of 11 levels of σ. In agreement with established theory, model cluster size PDFs near criticality (i.e., σ = 1) fit a −3/2 power law very closely (Fig. 2B, right, black) (Harris, 1989; Zapperi et al., 1995). Just as in the experiment, we computed κ based on CDFs of simulated spontaneous activity for different values of σ (Fig. 2B). We found that κ and σ were almost linearly related (Fig. 2D), which supports the following interpretation: In the experiments, κ ≅ 1 is close to criticality, κ < 1 identifies the subcritical regime, and κ > 1 is analog to the supercritical regime of the model.
Stimulus-evoked activity and dynamic range
After obtaining κ for a given experimental condition, we recorded the response R as a function of stimulus amplitude S (for peristimulus time histograms of response for different S, see supplemental material, available at www.jneurosci.org) Typical response curves from experiments and simulations are shown in Figure 3, A and B, respectively. We found that the shape of the response curves in the model closely matched the experimental findings. When excitatory synaptic transmission was reduced (κ < 1), the system was relatively insensitive (required a larger stimulus to evoke a given response). When inhibitory synaptic transmission was reduced (κ > 1), the system was hyperexcitable, with responses that saturate for relatively small stimuli. In the balanced E/I condition with κ ≅ 1, the range of stimuli resulting in nonzero and nonsaturated response was largest.
Maximal dynamic range at criticality, κ ≅ 1
For each response curve, we quantified the range of stimuli the network can process, i.e., the dynamic range Δ (see Materials and Methods). We found experimentally that Δ = 5.0 ± 0.1 (mean ± SE) under normal conditions, Δ = 2.4 ± 0.1 in the presence of PTX, and Δ = 3.4 ± 0.3 for AP5/DNQX (Fig. 3C) (F(2,44) = 11.3; p < 0.05 PTX and AP5/DNQX from normal). Importantly, the dynamic range was largest in unperturbed networks, in which neuronal avalanches are most likely to occur. These results were robust for different network sizes and different maximal stimulus amplitudes, whether or not the response curves reached saturation for all conditions (see supplemental material, available at www.jneurosci.org). Similar overall changes in Δ were also found in our simulations (Fig. 3D) (F(2,195) = 820; p < 0.05).
We then derived a tuning curve of Δ versus κ by combining all experimental conditions into one scatter plot (Fig. 4A). These data demonstrate that Δ is maximized and its variability is largest near κ ≅ 1. These findings agree well with our model including changes in Δ as the system is pushed away from κ ≅ 1, ∼10 dB drop (10-fold reduction in Smax/Smin) for a 30% change in κ (Fig. 4B). The tuning curve demonstrates that the change in the dynamic range of a network due to a shift in E/I depends on both the original, unperturbed state and the resulting change in κ.
Network tuning curve for dynamic range Δ near criticality. A, In experiments, Δ peaks close to κ ≅ 1 and drops rapidly with distance from criticality. Paired measurements share the same symbol shape; normal (no-drug) condition was measured just before the drug condition. Circles, Unpaired measurement. B, In simulations, Δ is also maximum for κ ≅ 1. Symbol indicates network size (circles, N = 250; squares, N = 500; triangles, N = 1000). Lines represent binned averages.
Discussion
We experimentally derived a tuning curve that linked the state of a cortical network with its ability to process stimuli. When the network was closest to criticality, as indicated by neuronal avalanches, κ was close to 1 and dynamic range was maximized. This is the first experimental demonstration to confirm theoretical predictions on the computational advantage of operating at criticality. Dynamic range has been predicted in simulations to peak at criticality (Kinouchi and Copelli, 2006). Our simulations advance previous studies by linking the dynamic range of a network with the spontaneous activity it generates. Because the dynamic range increases with the ability of a network to map input differences into distinguishable network outputs, our result is also closely related to network-mediated separation, which has been predicted to peak at criticality, at the transition from ordered to chaotic dynamics (Bertschinger and Natschlager, 2004; Legenstein and Maass, 2007). In contrast, our results show that variability of response to a given stimulus is highest at criticality. Further investigation of reliability versus variability in cortical networks is warranted.
Considering the simplicity of our model with all-to-all connectivity, absence of refractory period, and approximating inhibition by reducing σ, the overall agreement in the Δ–κ relationship between experiment and simulation is remarkable. The increase of variability in Δ as well as the drop in Δ due to deviation from κ = 1 was well matched between experiment and simulations. Such similarity supports the notion that universal principles, independent of system details, are found at criticality (Stanley, 1971; Bak and Paczuski, 1995; Jensen, 1998). The main quantitative difference was the lower Δ values for experiments compared to the model. Experimental noise, which is absent in the model, effectively adds a constant value to Smin and Smax, which systematically reduces Δ.
Further neurophysiological insight into our results can be gained from Figure 3. There it is shown that networks poorly discriminate small inputs in the hypoexcitable state, whereas they tend to saturate, failing to discriminate larger inputs in the hyperexcitable state. Both these reductions in performance result in reduced dynamic range compared to balanced networks. In line with these findings, dissociated cultures respond to inputs with a “network spike” if σ > 1 (Eytan and Marom, 2006) and display a “giant component” in a hyperexcitable regime, which reduces the ability to discriminate inputs (Breskin et al., 2006). The balance of excitation and inhibition has been shown to be crucial for developing cortical circuits to accurately process sensory inputs (Hensch, 2005). Our results suggest that, functionally, the balance of excitation and inhibition is achieved when the dynamic range is maximized and cortical networks operate at criticality.
Footnotes
-
This work was supported by the Intramural Research Program of the National Institute of Mental Health. R.R. acknowledges support from a Department of Defense Multidisciplinary University Research Initiative grant (ONR N000140710734) to the University of Maryland. T.P. is also grateful to the Swiss National Science Foundation (Grant PBEL2-110211) for financial support. We thank Craig Stewart for help with preparing the cultures.
- Correspondence should be addressed to Dr. Dietmar Plenz, Section on Critical Brain Dynamics, Laboratory of Systems Neuroscience, Porter Neuroscience Research Center, National Institute of Mental Health, Room 3A-100, 35 Convent Drive, Bethesda, MD 20892. plenzd{at}mail.nih.gov
