Abstract
Cognitive processes depend on synchronization and propagation of electrical activity within and between neuronal assemblies. In vivo measurements show that the size of individual assemblies depends on their function and varies considerably, but the timescale of assembly activation is in the range of 0.1–0.2 s and is primarily independent of assembly size. Here we use an in vitro experimental model of cortical assemblies to characterize the process underlying the timescale of synchronization, its relationship to the effective topology of connectivity within an assembly, and its impact on propagation of activity within and between assemblies. We show that the basic mode of assembly activation, “network spike,” is a threshold-governed, synchronized population event of 0.1–0.2 s duration and follows the logistics of neuronal recruitment in an effectively scale-free connected network. Accordingly, the sequence of neuronal activation within a network spike is nonrandom and hierarchical; a small subset of neurons is consistently recruited tens of milliseconds before others. Theory predicts that scale-free topology allows for synchronization time that does not increase markedly with network size; our experiments with networks of different densities support this prediction. The activity of early-to-fire neurons reliably forecasts an upcoming network spike and provides means for expedited propagation between assemblies. We demonstrate this capacity by observing the dynamics of two artificially coupled assemblies in vitro, using neuronal activity of one as a trigger for electrical stimulation of the other.
Introduction
Behaviors, from simple to most complex, involve orchestrated activation of neural cell assemblies. Functionally defined by Hebb, neuronal assembly is a group of cells that share similar static and dynamic response properties when activated through specific receptors, constituting “… the simplest instance of a representative process (image or idea)” (Hebb, 1949). In a series of classical electrophysiological studies (for review, see Mountcastle, 1998), as well as in later experiments in which large-scale imaging technologies were applied (Slovin et al., 2002; Ohki et al., 2005), the abstract notion of cell assembly was mapped to actual neural entities. These studies show that, depending on the nature and complexity of the stimulus, the numbers of neurons constituting an assembly range from hundreds to many hundreds of thousands (Roland, 2002; Derdikman et al., 2003). Neurons constituting an assembly synchronize during presentation of a matching stimulus, as well as during “ongoing activity,” in the absence of stimulation (Kenet et al., 2003). Regardless of the stimulus modality, stimulus complexity, or cortical area involved, the characteristic timescale of assembly activation is in the order of 0.1–0.2 s. This timescale emerges whether assembly activation is measured in the sensory (Super et al., 2001; Slovin et al., 2002), somatosensory (Derdikman et al., 2003), or motor (Riehle et al., 1997) areas. The same timescale characterizes the activation of “higher” cortical areas during categorization tasks (Keysers et al., 2001) and during pure internal events (Riehle et al., 1997). Thus, the 0.1–0.2 s timescale of assembly activation is a fundamental factor that might constrain the temporal aspects of cognition.
Because the timescale of assembly activation seems common to various cortical structures, it is conceivable that its biophysical characteristics may be obtained by experimentally analyzing a “generic” cortical neuronal assembly. A large set of experimental data and theoretical analyses indicates that a network of sparsely coupled cortical neurons developing in vitro might be useful as an experimental model for a generic assembly, keeping in mind the obvious constraints on extrapolations from in vitro to in vivo conditions (for review, see Corner et al., 2002; Marom and Shahaf, 2002). When developing on arrays of microelectrodes (Gross, 1979; Stenger and McKenna, 1994; Morin et al., 2005) (see Fig. 1a), the networks can be sampled simultaneously at many points and interrogated by site-specific stimuli with varying temporal and spatial structures, allowing a controlled biophysical examination.
In the present study, we use networks of cortical neurons developing in vitro on arrays of microelectrodes to describe the biophysical process underlying assembly activation in terms of population dynamics and topology. We show that assembly activation is a threshold-governed phenomenon. We denote the phenomenon “network spike” (NS). We point at the origin of the 0.1–0.2 s characteristic timescale and show that the underlying topology of effective connectivity is scale free. The impact of the dynamics and topology at the single assembly level on the synchronization and propagation of activity within and between assemblies is demonstrated experimentally.
Materials and Methods
Culturing the cortical neurons on multielectrode arrays.
Cortical neurons were obtained from newborn rats within 24 h after birth, following standard procedures (Shahaf and Marom, 2001; Marom and Shahaf, 2002; Eytan et al., 2003, 2004). The neurons were plated directly onto a substrate-integrated multielectrode array (MEA). The cultures were bathed in MEM supplemented with heat-inactivated horse serum (5%), glutamine (0.5 mm), glucose (20 mm), and gentamycin (10 μg/ml) and were maintained in an atmosphere of 37°C, 5% CO2/95% air in a tissue culture incubator as well as during the recording phases. Experiments were performed during the third week after plating, thus allowing functional and structural maturation of the neurons. MEAs of 60 Ti/Au/TiN electrodes, 30 μm in diameter, and spaced 200 μm from each other (Multi Channel Systems, Reutlingen, Germany) were used. The insulation layer (silicon nitride) was pretreated with poly-d-lysine. Experiments lasting >3 h were conducted using a slow perfusion system with perfusion rates of ∼100 μl/h.
Electrophysiological recordings.
A commercial 60-channel amplifier (B-MEA-1060; Multi Channel Systems) with frequency limits of 1–5000 Hz and a gain of 1024× was used. The B-MEA-1060 was connected to MCPPlus variable gain filter amplifiers (Alpha Omega, Nazareth, Israel) for additional amplification. Current stimulation through the MEA was performed using a dedicated eight-channel stimulus generator (Multi Channel Systems). The pair of electrodes for current passing were chosen based on practical considerations; in particular, we use electrodes that are not useful for recordings because of their noisy character. Stimulation parameters are detailed in the legends of Figures 6, 7, and 11. Data was digitized using two parallel 5200a/526 analog-to-digital boards (Microstar Laboratories, Bellevue, WA). Each channel was sampled at a frequency of 24,000 Hz and prepared for analysis using the AlphaMap interface (Alpha Omega). Thresholds (8× root mean square units; typically in the range of 10–20 μV) were defined separately for each of the recording channels before the beginning of the experiment. The data presented in the text is not spike sorted. The basic observation of nonrandom recruitment is completely preserved also after spike sorting using principal component analysis (supplemental Fig. 3, available at www.jneurosci.org as supplemental material). The supplemental figure also shows that each electrode in our setup senses approximately one to two neurons and rarely three neurons. The robustness of the main finding to spike-sorting procedure and the small number of neurons sensed by a single electrode led us to often interchange the terms “electrode activity” and “neuronal activity.”
Detecting network spikes.
Network spikes were detected using an activity threshold. A timestamp is defined when activity crosses a threshold of N number of action potentials recorded throughout the electrode array within a T millisecond time bin. In most cases (data of Figs. 1, 3a, 4b,c, 7, 9), N is equal to one-fourth of the (active) electrodes, and T = 3 ms. An active electrode is defined as demonstrating an average firing rate of >0.02 s−1 throughout the recording session. In several cases in which analyses required deviation from the above definition, threshold parameters are specified. The data around each timestamp (±300 ms) are extracted and stored. Temporal overlap between NSs was not allowed.
In our data, obtained by MEA with dimensions 1.96 mm2, no wave-like directional propagation is observed. This is consistent with data obtained by others that use that preparation (summarized by Marom and Shahaf, 2002) and is similar to prestructured preparations (slices and cultured slices) in which activity is not often wave like (Beggs and Plenz, 2004, their Fig. 2).
Results
Dynamics
Cortical networks in vitro are spontaneously active (Fig. 1b). The basic firing rate of individual neurons is 0.1–5 Hz, similar to figures obtained from cortical neurons in vivo. These values for spontaneous firing rates hold under continuous perfusion in vitro (see Materials and Methods); if perfusion is discontinued, the network gradually converges to a synchronic mode of action, with little (if any) activity between synchronizations. The spontaneous activity is completely abolished if excitatory synaptic transmission is blocked, indicative of the fact that spontaneous network activity is initiated by interactions at the population level (sporadic synaptic activity) rather than driven by self-pacing neurons (data not shown) (for review, see Marom and Shahaf, 2002). Over the past 20 years, physiologists repeatedly demonstrated that, like in vivo neurons (Kenet et al., 2003), cortical networks in vitro spontaneously synchronize once every 1–20 s (Fig. 1b,c), generating assembly activity events (Habets et al., 1987; Ramakers et al., 1990; Corner and Ramakers, 1991, 1992; Muramoto et al., 1993; Maeda et al., 1995, 1998; Kamioka et al., 1996; Nakanishi and Kukita, 1998; Ben-Ari, 2001; Corner et al., 2002; Marom and Shahaf, 2002; Beggs and Plenz, 2003, 2004; Wagenaar et al., 2005, 2006). Synchronizations with similar kinetics may also be evoked by site-specific electrical stimuli (Jimbo et al., 1999; Eytan et al., 2003). These synchronous activities, represented in terms of the number of action potentials recorded in millisecond time bins, are here referred to as network spikes (Fig. 1c,d). Results obtained by systematically varying the threshold for network spike detection indicate that, in most cases, once a network spike starts to assemble, the totality of active electrodes in the preparation is eventually recruited. In this respect, and based on the nature of responses to external stimulation (demonstrated later in Fig. 6c, inset) (Jimbo et al., 1999; Eytan et al., 2003), network spikes have the flavor of “all-or-none” phenomena (Crain, 1976). Figure 2, a and b, demonstrates this point, showing the number of different electrodes participating in network spikes (color coded) as a function of the threshold for NS detection (the two panels summarize results from two different networks). Note, however, that a closer examination of these histograms (Fig. 2b in particular) reveals the existence of “aborted” network spikes, most readily detected using reduced thresholds. To further clarify the picture, ∼900 network spikes that served for the construction of Figure 2b are sorted and plotted in Figure 2c: each column depicts the identity of the electrodes (y-axis) that participated in a given network spike (x-axis). The events are ordered such that the group of aborted network spikes is clearly seen. Note that the electrodes recruited by aborted spikes constitute a subpopulation of those recruited by fully developed spikes; we will return to this point later.
a, Cortical network on substrate-embedded multielectrode array. The dark circle is a 30-μm-diameter electrode. Neurons are tagged using green fluorescent protein. b, Example of spontaneous activity simultaneously recorded from eight different channels. Top, At 500 s. Bottom, Higher temporal resolution of 30 s from the top panel (extracted section is depicted by a dark bar). A box marks a single event of synchronous activity. c, Top three traces show examples of individual synchronous events in terms of number of spikes recorded in 60 electrodes (1 ms time bins). The average of 273 such events (NSs) is shown. d, Example of average NSs recorded over 1 h from different networks (normalized amplitudes).
a, b, Color-coded representations of the number of different electrodes participating in a network spike (x-axis) as a function of threshold for detection of network spikes (y-axis). Colored scale bar (on the right) depicts number of occurrences. The two panels show distributions obtained over a 2 h period in two different networks that faithfully represent the behavior of the entire set of ∼40 networks served for this study. Threshold is expressed in terms of the number of electrodes that are required to be active (within a 3 ms time bin) for time-stamping an NS. Time bin width was also subjected to systematic variations (from 3 to 30 ms) with no qualitative effect (data not shown). c, Approximately 900 NSs obtained over 1 h of recording from the network shown in b, using a very low threshold (4 action potentials in 10 ms time bin). The identity of the electrodes participating in each NS is shown (black). Although in the majority of the NSs all active electrodes appear, there exists a subpopulation of aborted NSs in which only a subset of the active electrodes participate. Note that the NSs are not depicted according to the chronology of their appearance; rather, for purposes of clarity, they are sorted using a hierarchical clustering algorithm.
Spatial sampling is limited in our experimental setup to 60 substrate-embedded microelectrodes; however, averaging enhances analysis of the kinetics underlying a network spike. Averaging thousands of spontaneous network spikes from 21 experiments reveals that the early phase of the spike, in which the activity of the assembly just starts to increase, fits an exponential population growth model (Fig. 3a). The rate of recruitment, σ, defined as the number of active electrodes at a given point in time A(t), divided by the number of active electrodes at (t − Δt), is 1.045 (range of 1.02–1.07; SD of 0.015; Δt = 1 ms; 21 networks, 5796 NSs). Thus, the pool of active electrodes increases by ∼5% with each millisecond. The rate of recruitment seems to be independent of the network spike size; a Pearson’s correlation coefficient of r = −0.14 was calculated between σ and spike size (ranging from 7 to 17 active electrodes/ms; 21 experiments). Furthermore, Figure 3b (top) shows a distribution of network spike sizes for one given network, obtained using a detection threshold of four action potentials per 10 ms time bin; recruitment rates calculated separately from the right and left sides of the distribution are very similar (1.05 and 1.04, respectively). The slight tendency for faster recruitment in larger network spike sizes, manifested by both the negative Pearson’s coefficient and the exponential functions of Figure 3b (bottom), is addressed later.
a, Inset, Grand average network spike, calculated by averaging 21 average NSs similar to those shown in Figure 1d (total of 5796 NSs). Marked in brown is the initial segment, enlarged (brown dots) in the main figure. An exponential growth equation, A(t) = a + b × e(σ−1)t, was fitted to this initial segment; the resulting function (a = 0.05, b = 0.01, σ = 1.045) is depicted by a continuous black line. b, Top, Distribution of NS peak activity for one network that demonstrated a relatively broad spectrum of NS amplitudes. A very low threshold for NS detection was used here (4 action potentials in 10 ms bin). Bottom, Early recruitment phases for 207 low-amplitude NSs (≤6 action potentials/ms; depicted in black) and 208 high-amplitude NSs (≥10 action potentials/ms; depicted in gray), fitted with exponential growth equations, yielding σ values of 1.04 and 1.05, respectively.
We estimated the rates governing the kinetics of network spikes by comparing recordings made in the presence and absence of bicuculline, a pharmacological agent that blocks fast inhibitory synaptic transmission. Figure 4a shows that, in the presence of 5 μm bicuculline, the number of aborted network spikes dramatically decreases, and the activity after the peak of an individual network spike becomes more vigorous. In Figure 4b, an average network spike in the presence of bicuculline (depicted in brown) is superposed on top of the control average network spike, showing that the inhibitory subnetwork is predominantly affecting the late phase of relaxation from an NS and has a characteristic timescale of ∼30 ms (Fig. 4b, bottom right inset). Indeed, in the absence of an inhibitory effect (brown trace), the network spike does not relax completely; the residual activity is sustained for several hundred milliseconds. Figure 4b also shows that the initial decline of activity from the peak of an NS to baseline is dominated by cellular-level processes: much of the force that restores the level of activity after the peak of an NS is insensitive to bicuculline, suggestive of the role played by synaptic depression, refractoriness, and cellular adaptation in this phase of the NS. Analysis of the impact of bicuculline on the rate of recruitment, σ, was enhanced by collecting recruitment phases before and after addition of 5 μm bicuculline to the bath medium. This procedure, repeated in 16 different networks and summarized in Figure 4c, reveals an acceleration in the rate of recruitment when the inhibitory subnetwork is pharmacologically neutralized, leading to σ values in the range of 1.2/ms.
a, Network spikes recorded over a period of 2 h in the absence (left) and presence (right) of 5 μm bicuculline. These NSs were obtained using a threshold of three action potentials within a 10 ms time bin. Each black dot marks an action potential detected in any of the electrodes during the ±250 ms surrounding the NS threshold. NSs are ordered using a clustering algorithm to enhance visualization of bicuculline effects: a reduced number of aborted NSs and a more vigorous activity within each NS. b, An average NS obtained in the presence of 5 μm bicuculline (brown). Average NS in control solution is shown in black. Top right inset, Fitted fast declining phase in the presence of bicuculline, representing the timescale of cellular-level (intrinsic) restoring forces. Bottom right inset, Fitted slow declining phase of the control NS, representing the timescale of restoring force acting through the inhibitory subnetwork. Note the timescale separation between the two types of restoring forces involved (6.5 ms for the effect of cellular-level forces; 33.7 ms for the effect of the inhibitory subnetwork). c, Averaged phases of recruitments extracted from 16 experiments before (black) and after (brown) the addition of 5 μm bicuculline to the bathing solution. Each recording episode lasted 1–2 h. Broken lines depict 1 SD (one-sided for clarity). Thick lines depict segments for which a single-exponential recruitment function could be reliably fitted for both conditions; σ values obtained were 1.06 and 1.22 for control and bicuculline conditions, respectively.
The kinetics of recovery from the effect of restoring forces (“refractory period” of an NS) occurs on a seconds timescale and is not markedly affected by blockade of inhibitory synapses. These recovery kinetics were estimated by observing indicators for network excitability after a network spike. In particular, we used electrodes that fire steadily at a relatively high frequency between network spikes (Fig. 5); immediately after a network spike, they become silent (Fig. 5a,b) as a result of the effects of cellular- and network-level restoring forces. Thus, the kinetics of recovery to the “ongoing” firing rate in these electrodes, after a network spike (Fig. 5b, arrows), represents the kinetics of recovery from the effects of restoring force. Obtaining average recovery times in the presence and absence of bicuculline (Fig. 5c) provides an estimate of recovery kinetics of intrinsic (cellular adaptation, refractory period, and synaptic depression) and extrinsic (inhibition) restoring forces. We find that, under both conditions (with and without bicuculline), the recovery timescale is within the range of several seconds (Corner and Crain, 1972): the recovery timescale in the control experiment is 5 s (average of 96 recovery events). In three bicuculline experiments, the obtained recovery timescales are 2, 2.6, and 8.6 s (1069, 828, and 726 recovery events, respectively).
a, b, Firing rate of a relatively active neuron; the rate increases during an NS, decreases dramatically after the NS, and then gradually recovers (arrows in b). These kinetics are used for estimation of recovery from a network spike. c, Fitted recoveries, such as those depicted by arrows in b, in control (black) and bicuculline (brown) solutions. [Note that the inter-NS intervals in bicuculline solution are short; points at which the average firing rate during recovery was “contaminated” by the presence of NSs were omitted (bins 7–9, 11, 14 s)].
The above kinetic analyses provide much of what is required to heuristically model the network spike phenomenon in terms of standard population dynamics. We start by pointing out that the 0.1–0.2 s timescale of assembly activation reemerges when the recruitment rate, σ = 1.04 ms−1, is embedded into a normalized logistic growth equation, dA/dt = (σ − 1) A (1 − A), that describes the fraction of active electrodes (A) as a function of time. The term (1 − A) constrains the model to a finite population size. The logistic equation is widely used in population models in which the rate of growth is limited by the carrying capacity (i.e., maximum size) of the population; thus, it can only account for the rising phase of the network spike. To account for the entire shape of a network spike, including its relaxation, the impact of restoring forces at both the cellular level (e.g., adaptation, refractoriness, and synaptic depression) and the network level (activation of inhibitory subnetwork) need be considered. To that end, we describe the population growth by a modified version of the logistic equation dA/dt = (s(A,t) σ − 1) i(A,t) A (1 − A). s(A,t) and i(A,t) are kinetic variables ranging from 0 to 1, embodying cellular and network-level restoring forces, respectively. The assumed kinetic schemes for s(A,t) and i(A,t) are similar to the form used by Hodgkin and Huxley (1952) in their description of a neuronal action potential; that is, (1 − i) ⇔ i and (1 − s) ⇔ s, with transition rates being exponential functions of A(t). The underlying assumption is that the recruitment of inhibition and the processes of cellular and synaptic adaptation are nonlinearly dependent on activity. To avoid discontinuities in the model, we used exponential terms; attempts to physically interpret these exponential terms will be presumptuous, at present. Figure 6 summarizes our kinetic analyses in the form of a numeric solution. In Figure 6a, the time constants used for simulating i(A) and s(A) are shown. The rate equations for i and s were devised to approximately match the timescales of activation and recovery of the restoring forces at the extremes of A = 1 and A = 0, based on data shown in Figures 3⇑–5. Figure 6, b and c, shows spontaneous and evoked network spikes, respectively, generated by the model with Gaussian noise added. The inset to Figure 6c shows data from an experiment in which a network was excited by applying short biphasic current pulses between a pair of electrodes. Taken as a whole, the modified continuous version of the logistic equation faithfully reconstructs key features of the experimental records, yielding a network spike with a 0.1–0.2 s characteristic timescale, which is fairly insensitive to inhibition in its present form.
Numerical solution of the modified logistic equation dA/dt = (s(A,t) σ − 1) i(A,t) A (1 − A) for spontaneous and evoked NSs. a, Time constants of i(A) and s(A). Exponential rate equations for i and s were devised to approximately match the timescales of activation and recovery of the restoring forces at the extremes of A = 1 and A = 0, based on data shown in Figures 3⇑–5. The equations used the following (time in millisecond units): βs = 30e−20(A(t) + 0.5); αs = 0.00001e10A(t); βi = 10e−20(A(t) + 0.5); αi = 0.000001e10A(t). 0 ≤ A(t) ≤ 1; σ was set to 1 below a threshold of A(t) = 0.05. Once the threshold is crossed, σ is set to a value ranging from 1.04 to 1.06. Noise was generated using normal (Gaussian) distribution around mean A(t), with SD of 0.001 and 0.0002 for b and c, respectively. For the case of evoked network spikes (c), a fraction of the population is excited (i.e., stimulated) from around threshold and higher in discrete steps. Inset in c shows data from an experiment in which a network was excited by applying short (0.4 ms) biphasic current pulses between a pair of electrodes [stepping from 10 μA (bottom trace) to 80 μA (top trace)]. Average evoked responses to 15 presentations of stimuli at each amplitude are shown (each trace is 300 ms long). Missed responses were not included except for the case of 10 μA stimulation amplitude. We observed 15 of 15 missed responses for 10 μA stimuli, 5 of 15 for 20 μA, 4 of 15 for 30 μA, and 2 of 15 for 40 μA; beyond 40 μA stimulation, no missed responses were observed.
Effective topology
Measurements in vitro (Ikegaya et al., 2004; van Pelt et al., 2004) as well as in vivo (Abeles, 1991; Tsodyks et al., 1999; Buzsaki et al., 2004) indicate that the sequence of neuronal activation within a synchronization event is nonrandom and strongly constrained by the pattern of population activity. Figure 7 shows that, indeed, the recruitment of neuronal activity within the NS is hierarchically structured; there exists a stable subset of privileged neurons that reliably increase their firing rates tens of milliseconds before the peak of the network spike. The two plots of Figure 7a summarize data collected over 1 h each (averages of >200 NSs in each). Note that the pattern of activity obtained by averaging NSs during the first hour of recording is very similar to the pattern obtained over the seventh hour; such stability is in agreement with data published by van Pelt et al. (2004). Particularly relevant for the subject matter of the present study is the fact that the same privileged neurons appear in both recording sessions. As shown in Figure 7b, the privileged neurons within a given assembly predict an upcoming network spike regardless of its activation source: whether the network spike starts spontaneously (left panel) or is evoked externally by applied electrical stimuli at different sites (right panel), the same privileged neurons are the first ones to be activated. Furthermore, as suggested by the results of Figure 8, the same privileged neurons are the first ones to be active whether the synchrony evolves to a fully developed network spike or an aborted one. The activity of privileged neurons reliably predicts an upcoming network spike as early as 100 ms before the peak of the spike. Figure 9a demonstrates this point by showing the cumulative number of spikes emitted by six neurons in a given network, some of which are privileged.
a, Firing probability of neurons as a function of time surrounding an NS. Color scale ranges from 0 to 1 spike per 5 ms (the few cases in which >1 spike occurred during a 5 ms bin are represented as 1 spike/bin). b, In a given network, the early-to-fire neurons are similar for spontaneous (left) and evoked (right; 3 different stimulation sources) NSs. For the case of spontaneous NS, arrows point to times for which firing probabilities are presented. Probabilities of firing 0–5 and 25–75 ms after stimulation are shown for three different stimulation sites (S1, S2, S3). Short (0.4 ms) biphasic 30 μA current pulses between each of three different pairs of electrodes were applied. Average responses to 60 presentations of stimuli for each stimulation point are shown. Horizontal lines depict four examples of early-to-fire neurons. Poststimuli time histograms for the three stimulation sites are shown at the bottom right.
A network with a relatively large number of aborted network spikes is subjected to analyses aimed at comparing early-to-fire neurons in aborted versus full-blown NSs. a, All 1087 NSs detected over 1 h of recording are shown in a raster plot, ± 250 ms surrounding a detection threshold of four action potential in 10 ms time bin. A black dot depicts an action potential recorded in any of the electrodes during the NS. Responses were reordered using a clustering algorithm to enhance visual separation between full-blown and aborted network spikes. b, The distribution of peak number of action potentials in each NS (1 ms time bin; same procedure as in Fig. 3b). Distributions around threshold for 537 NSs with a peak of ≤5 action potentials (c) and 437 NSs with a peak of ≥10 action potentials (d) are shown using color scale ranges from 0 to 0.33 action potentials per 5 ms (the cases in which >0.33 action potentials occurred during a 5 ms bin are represented as 0.33 action potentials per bin). Note that the same early-to-fire neurons are active for both subsets of NSs. The correlation between firing rates of the electrodes of both subsets at times ranging from −300 to −50 ms before NS threshold is 0.98.
a, Cumulative number of spikes emitted by six neurons in a network, as a function of time relative to the occurrence of an NS (for definition of NS timestamp, see Materials and Methods). Inset, Distribution of intervals between spontaneously occurring NSs (491 synchronous events recorded from this network). b, Firing rate distribution 100–75 ms before the NS (black) and throughout the NS (gray; 1200 neurons from 20 networks; 9400 NSs). The distribution is fitted by a power-law function with a scaling power d ≈ −2. Inset, A total of 387 different neurons (12 different networks) that were active as early as −100 to −50 ms before an NS threshold crossing point were selected. For these neurons, the coefficient of variance (CV) of instantaneous frequency (manifested by interspike interval) during 50 ms after threshold crossing was calculated and plotted as a function of their firing rate ±300 ms surrounding the NS detection threshold. Note that the coefficient of variance is negatively correlated to the neuronal firing rate (r = −0.22; p < 0.0001).
Recall that none of the neurons in the network fires unless it is driven by other neurons. This means that the firing rate of an individual neuron in general, and during the recruitment phase in particular, reflects its sensitivity to the activity of other neurons in the network. Note that this sensitivity is mediated by a combination of factors, including, for instance, the number of synaptic inputs that the neuron receives, the distribution of synaptic weights in its dendritic tree, resting potential, firing threshold, dendritic conductances, and dynamics of cellular adaptation. Thus we use the term effective connectivity (rather than simply “connectivity”) to designate the sensitivity of a neuron to network activity as reflected in its firing rate. Histograms of firing rates obtained from 1200 active electrodes (20 networks) throughout 600 ms that surround an NS, as well as within the time window of −100 to −75 ms before its peak are shown in Figure 9b (gray and black circles, respectively). Both distributions are broad, and the latter is fitted by a power law over ∼2 orders of firing rate magnitudes (covering the physiological range of firing rates). The power-law distribution shown in Figure 9 does not artificially result from the variance in average firing rates of the networks from which neuronal activities were pooled. To establish this point, we looked at firing rate distributions obtained from two subgroups of networks. In one subgroup (six networks), the average electrode firing rate was 9.9–13.7 Hz during a 600 ms time window surrounding the network spike. In the other subgroup (seven networks), the average firing rate was lower (3.4–5.2 Hz). One-way ANOVA suggests that the networks, within each of the two subgroups, do not significantly differ in the firing rate distributions (all 60 electrodes from each network were taken into account, including those that did not show any activity; F ratio = 1.25, p = 0.28 and F ratio = 1.01, p = 0.41, for the high- and low-firing-rate networks, respectively). Supplemental Figure 1 (available at www.jneurosci.org as supplemental material) shows that the firing rate distributions of electrodes from both low- and high-firing-rate subgroups are broad and described by a power-law function (d = −1.85, R2 = 0.87 and d = −2.01, R2 = 0.91, respectively). Analysis of firing rate fluctuations (Fig. 9b, inset) indicates that the coefficient of variance for strongly active electrodes is reduced compared with active electrodes that fire less during the recruitment phase. Assuming that the dynamics of individual synapses as well as the membrane properties of neurons that fire more are not different from those that fire less, the inset of Figure 9b implies that neurons that fire more integrate over more synaptic input compared with neurons that fire less.
Broadly distributed connectivity in general, and power-law distributed connectivity in particular, are usually associated with random graphs in which the average topological distance between nodes increases very slowly with the number of nodes, despite a large local interconnectedness (Strogatz, 2001). In accordance, under certain conditions that seem applicable to neuronal networks (Nishikawa et al., 2003; Motter et al., 2005), compared with Erdös-Rényi type of connectivity, dynamic systems coupled in this way display enhanced signal propagation speed that increases very slowly (Strogatz, 2001; Newman, 2003) or even decreases (Barthelemy et al., 2004) as the network size gets bigger. To test this prediction, we prepared networks at two different neuronal densities (0.1× and 2× the standard density). The difference between the densities of the two groups as they matured (third week in vitro) did not remain as large as was when the networks were plated, probably because of compensatory mechanisms that are beyond the scope of the present study. However, there was still a very clear difference between the neuronal densities of the two groups, observed structurally (Fig. 10a) and functionally (Fig. 10b). Figure 10c shows that indeed the time it takes for the network to synchronize does not increase as the number of active electrodes participating in the synchrony increases; if anything, the denser the network is, more active electrodes participate in the synchrony and synchronization becomes faster.
Effect of neuronal densities on synchronization time. a, Photographs of exemplar mature networks (3rd week in vitro) plated at 0.1× (left) and 2× (right) the standard density. Exemplars of the resulting network spike amplitudes are shown in b and normalized in the inset in c, demonstrating that the time it takes for the network to synchronize does not increase as the number of neurons participating in the synchrony increases. To allow analyses under conditions of low density, a 10 ms time bin for threshold detection was used for both (high and low) densities. c, A summary of 13 different networks, showing that the time-to-peak of the average network spike decreases as the network becomes more dense. Time-to-peak was calculated from 10% activity to peak activity.
Hierarchical recruitment within an assembly provides a powerful mechanism for modulation of time delays between coupled assemblies. To demonstrate this feature, we electrically coupled a pair of assemblies, X → Y, using a stimulus generator. Obviously, under such artificial conditions the “actual” time delays are arbitrarily dictated by setting the coupling stimulus parameters. It would be expected that “ranking” the capacity of different neurons from assembly X to speed up the synchronization between X and Y is independent of the artificial stimulation parameters. Figure 11, a and b, shows that, by clamping the amplitude and locus of the input stimulus to Y, while changing the identity of the triggering neuron (in X), the time delays between network spikes in X and Y may be strongly affected; different neurons, by virtue of their effective connectivity, may consistently yield different time delays. When assembly Y “reads” the activity of assembly X through poorly connected neurons, the time delay between the network spikes in X and Y is large. However, when highly connected neurons are read, the network spike in Y can appear simultaneously with that of X or even precede it.
a, NSs of artificially coupled assemblies. Activity of different neurons in assembly X (brown) is used to trigger the delivery of a fixed current amplitude stimulus (0.4 ms, biphasic 50 μA) to Y assembly (black). Different traces represent the impact of using different triggering neurons from X. b, Summary of results from four artificially coupled pairs of networks (4 different symbols). Time delays between average NSs in X and Y are shown as a function of the triggering neuron effective connectivity. As an index for the effective connectivity in this figure, we used the average number of action potentials emitted by a neuron over the time period from −300 to −15 ms relative (i.e., before) to an NS timestamp. Note that values <1 indicate that a triggering neuron does not emit action potentials before each network spike during the designated period, in which case, it is its participation in the NS at later stages that causes assembly Y to be ignited.
Discussion
Using a substrate-integrated multielectrode array, we measured and characterized the kinetics of neuronal assembly activation in vitro. Expressed in terms of population firing rate, assembly activation is an all-or-none-like threshold-governed phenomenon that is very reminiscent of an action potential in a single neuron, hence the term network spike. The recruitment of neurons during the early phases of a network spike follows a single-rate exponential process that amounts to a characteristic timescale of 0.1–0.2 s for achieving synchronization. The measured recruitment rate (σ) is close to unity. Beggs and Plenz (2003, 2004) used a related measure to describe the propagation of local field potentials in cortical slices. Basing their elegant analyses on equations that govern avalanches, they also report a branching parameter close to unity, a value that optimizes information transmission in feedforward networks but prevents runaway network excitation. Although the fit of a single-exponential recruitment to the average data does seem tight (Fig. 3), small variations around the value of σ (SD of 0.015) are not negligible, as one might be tempted to assume. Incorporating σ values ranging from 1.025 to 1.055 (±1 SD around σ = 1.04 ms−1) into a logistic growth equation yields a range of ∼80 to ∼180 ms in the time to half-maximum activation, respectively. Indeed, examination of the traces shown in Figure 1d suggest that such a range is realistic.
The all-or-none nature of network spikes reported here seems incongruent with Beggs and Plenz (2003) power-law size distribution of avalanches (i.e., bursts of activities). Moreover, analysis of our data using the procedure of Beggs and Plenz results in bimodal distributions for both avalanche size (number of participating electrodes) and avalanche duration (supplemental Fig. 2, available at www.jneurosci.org as supplemental material). Our bimodal distributions are very similar to the bimodal avalanche size distributions seen in the presence of picrotoxin, a fast inhibition blocker (Plenz, 2005). We believe that the difference between the two observations is attributable to the fact that Beggs and Plenz were using a preparation cut out from a prestructured network with conserved layer organization; we, conversely, use an ex vivo spontaneously developing preparation. Indeed, a recent abstract of Plenz’s group (Stewart et al., 2005) indicates that the avalanches they observe originate in superficial layers and are also confined to these layers most of the time.
The relative impact of cellular- and network-level restoring forces was estimated using pharmacological manipulations. Although later phases of the network spike are affected to some extent by activation of the inhibitory subnetwork, the major restoring force that operates throughout a network spike arises from cellular processes [most probably synaptic depression (Eytan et al., 2003)]. A logistic equation with kinetic variables that represent the various restoring forces accounts for the 0.1–0.2 s time-amplitude trajectory of a network spike and is offered as a framework for mathematical modeling of assembly synchronization events.
The contribution of different neurons to the various phases of the network spike is not random. The phase of the network spike within which a neuron fires is strongly constrained; some neurons consistently fire at the very early phase of the network spike, whereas others start to fire at later phases of the spike. Nonrandom, sequenced activation of neurons that is strongly constrained by the pattern of population activity was demonstrated in vitro (Ikegaya et al., 2004; van Pelt et al., 2004), as well as in vivo (Abeles, 1991; Tsodyks et al., 1999; Buzsaki et al., 2004). Grinstein and Linsker (2005) addressed the role played by network topology in determining the nature of synchronous neural activity. They show that, in contrast to random Erdös-Rényi networks, networks having a power-law connectivity distribution generate large synchronous firing peaks dominated by a small subset of nodes. Here we used the firing rates of individual neurons during the early recruitment phase as indicators for their effective connectivity. We justify our approach by having observed that neurons in our preparation fire only in response to synaptic input, not spontaneously; hence, firing that starts earlier means more sensitivity to network activity, that is, more effective convergence. Firing that starts earlier also probably means more effective divergent connectivity in the sense that the more a neuron fires, the higher its impact on the propagation of activity and assembly synchronizability. We show that the distribution of firing rates during the early phases of the network spike is broad and described by a power law, which would be consistent with an essentially scale-free topology of connectivity.
It is tempting to think of the above results in the context of observations in vivo. The time to reach synchronization in our in vitro preparation (that contains ∼100,000 neurons) is 0.1–0.2 s. This is also the time it takes for a single cortical column to reach synchronization (Derdikman et al., 2003), as well as for whole brain areas such as primary visual area V1 and secondary visual area V2 (Slovin et al., 2002), and in fact for the entire brain; indeed, 0.1–0.2 s is also the accepted figure for simple behavioral reaction times. To a large extent, it seems that the time to reach neural synchrony is size invariant. In recent years, the relationships between network size and time to synchronize is intensively studied as part of the growing interest in the field of complex random graphs (for review, see Newman, 2003). A wide range of biological and social networks (graphs) has very broad connectivity distributions; these graphs are often referred to as “scale-free” networks, designating the power-law distribution of connectivity that is characteristic of some of them. In such cases, the average topological distance between nodes increases slowly with the network size, despite a large local interconnectedness (Strogatz, 2001). In the neurobiological context, scale-free topology within neuronal networks was suggested to provide an economical means for global synchrony and oscillations at multiple timescales (Buzsaki et al., 2004; Sporns et al., 2004; Grinstein and Linsker, 2005). As explained above, under certain conditions that might be applicable to neuronal networks (Nishikawa et al., 2003; Motter et al., 2005), dynamic systems coupled in this way display enhanced signal propagation speed that is invariant or even decreases with network size (Strogatz, 2001; Newman, 2003; Barthelemy et al., 2004). We suggest thinking in such terms on the mechanism of size invariance in the time it takes to synchronize networks in the brain. Our contribution to this line of thought is to show that increasing the number of neurons that participate in a network spike causes a decrease in time to reach synchrony.
Hierarchical recruitment of neurons during synchrony provides a means for modulation of time delays between sequential activations of coupled assemblies. The reliability of forecasting a network spike tens of milliseconds before its peak based on the activity of early-to-fire neurons, together with the slow integration toward a network spike, allows a wide temporal range for such modulation. Assembly coupling, presented in Figure 11, although limited by the artificiality of the stimulus, demonstrates such temporal modulation: the time delays between the activities of two coupled assemblies may be dramatically reduced. This scheme provides a possible explanation for the enigmatic rapidity of processing in pipelines of neuronal assemblies. Consider, for instance, the case of visual categorization experiments. Electrophysiological and psychophysical analyses show that the visual system recognizes a category in less than 
of a second (Keysers et al., 2001; VanRullen and Koch, 2003). Taking into account the number of neural cell assemblies through which the signal travels, such rapid processing implies 0.01 s of activity within each assembly before it is forwarded to the next one. This number is difficult to reconcile with the 10 times slower process (0.1–0.2 s) required for a full-blown activation of an individual assembly in response to a stimulus. Figure 11 shows that, in principle, when relying on the activity of early-to-fire neurons in the hierarchy of recruitment within an assembly, such short response times are accounted for (Thorpe and Fabre-Thorpe, 2001; Thorpe, 2002). Figure 11 also suggests that, in principle, a driven (next-in-line) assembly may synchronize simultaneously with, or even before, its driving (forerunner) assembly. This possibility surfaces a potential methodological difficulty in the interpretations of functional macroscopic neural data (e.g., electroencephalogram, magnetoencephalography, and functional magnetic resonance imaging) that rely on time delays between activation of assemblies. For instance, when two assemblies (X and Y) seem to be active together, the possibility of X driving Y or vice versa would be ruled out; it goes without saying that, when macroscopic activity from assembly X precedes that of Y, concluding that the latter being a cause of the former would seem unreasonable. The results of Figure 11 imply that, in constructing schemes of activation paths, much care should be exercised when relying on macroscopic time delays.
Footnotes
-
This work was partially supported by grants from the Israel Science Foundation, the National Institute for Psychobiology in Israel, and the Minerva Foundation. We thank Ella and Vladimir Lyakhov for technical assistance and Dr. Noam Ziv and Amir Minerbi for preparation of photo images.
- Correspondence should be addressed to Dr. Shimon Marom, Department of Physiology and Biophysics, Faculty of Medicine, Technion–Israel Institute of Technology, P.O. Box 9697, Haifa 31096, Israel. Email: marom{at}tx.technion.ac.il
