## Abstract

The synfire chain model has been proposed as the substrate that underlies computational processes in the brain and has received extensive theoretical study. In this model cortical tissue is composed of a superposition of feedforward subnetworks (chains) each capable of transmitting packets of synchronized spikes with high reliability. Computations are then carried out by interactions of these chains. Experimental evidence for synfire chains has so far been limited to inference from detection of a few repeating spatiotemporal neuronal firing patterns in multiple single-unit recordings. Demonstration that such patterns actually come from synfire activity would require finding a meta organization among many detected patterns, as yet an untried approach. In contrast we present here a new method that directly visualizes the repetitive occurrence of synfire activity even in very large data sets of multiple single-unit recordings. We achieve reliability and sensitivity by appropriately averaging over neuron space (identities) and time. We test the method with data from a large-scale balanced recurrent network simulation containing 50 randomly activated synfire chains. The sensitivity is high enough to detect synfire chain activity in simultaneous single-unit recordings of 100 to 200 neurons from such data, enabling application to experimental data in the near future.

## INTRODUCTION

Although a tremendous amount of detail exists about properties and connectivity of cortical neurons (Braitenberg and Schüz 1991; Douglas and Martin 2004; Petersen and Sakmann 2000; Song et al. 2005; Yoshimura et al. 2005) the neuronal organization that underlies computational processes is poorly understood. Neurons usually do not act alone and, since Hebb (1949), the idea of a neuronal assembly has become widely accepted as a functional unit. There is, however, a continuing debate in the literature over whether representation of information in and transmission of information between such assemblies is in the form of firing rates or in the form of spatiotemporal firing patterns. One particularly intriguing suggestion for assembly organization is the synfire chain first proposed by Abeles (1991). The original version was a chain of length *l* composed of groups of *w* neurons each (the synfire chain width), with full unidirectional connectivity between successive groups. A group together with its output connections can interchangeably be called a link in the synfire chain. Such a chain structure propagates near-synchronous activity in each group along successive groups like a row of dominoes. Detailed properties of various versions of such systems have received considerable theoretical attention; stability properties (Diesmann et al. 1999) have been studied as well as structural variations such as partial group to group connectivity or such as feedback when a given neuron occurs in more than one group. There also have been studies of partially interconnected synfire chains to define the conditions in which activity can jump between chains (Hayon et al. 2004, 2005). Recent theoretical studies suggest that coupled synfire chain structures can demonstrate compositionality, the hierarchical representation of complex entities in terms of parts and their relations (Bienenstock 1995; Hayon et al. 2004; Schrader et al. 2007).

Intriguing though the synfire concept may be, to date there has been no direct way to examine whether such structures occur in real brains. One approach might be to look for the repeating multineuron firing patterns that would be the signature of synfire activity. There are good methods (Abeles and Gerstein 1988; Gerstein 2004) to detect and test for significance any repeating spatiotemporal firing patterns that exceed the number expected by chance. Such patterns have been found in repeating behavioral task situations (Baker and Lemon 2000; Baker et al. 2001; Kohn and Smith 2005; Prut et al. 1998) and, in several papers, claimed as evidence for synfire chain activity (Abeles et al. 1993; Ikegaya et al. 2004). With parallel recording of the order of 10 neurons the likelihood that the small sample includes several neurons either from the same or successive links of a particular chain is vanishingly small, so that such patterns are more likely to directly represent information coding and flow.

However, as massively parallel recordings become available (electrode methods: Csicsvari et al. 2003; Donoghue 2002; Euston et al. 2007; Nicolelis and Ribiera 2002; Rolston et al. 2007; Warren et al. 2001; optical methods: Ohki et al. 2005, 2006; Sasaki et al. 2007) the sampling problem is eased. The pattern program (Abeles and Gerstein 1988) can easily deal with the larger number of neurons and will detect millions of candidate repeating patterns, some of which will be significant by comparison to surrogate data (proper choice of surrogate is a nontrivial problem in its own right). At this stage synfire chains can only be inferred (Abeles et al. 1993; Ikegaya et al. 2004), since they are a possible but not unique source for patterns. It is then necessary to sort, collate, and variously compare the large number of significant detected patterns to gradually build up a meta organization like several synfire chains. The process would be like assembling the pieces of a jigsaw puzzle, in this case making appropriate matches and interleavings of pattern occurrence time, neuron identities, and spike time interval structure. Such programs have not yet been tried and will have unknown and probably data dependent computing time requirements.

In this study we develop a new method that can rapidly and clearly detect synfire chain activity and identify the participating neurons in a massive recording of even thousands of neurons, but that also is sensitive enough so that it continues to detect and identify even when ≤200 neurons are available. We demonstrate the method's comparative run-time requirements and its effectiveness for detection and for eliciting some aspects of chain connectivity using data from a large-scale cortical network model with and without embedded synfire chains.

## METHODS

### Simulation

Development and evaluation of a new tool for spike train analysis can only be done using data with well-defined properties. The data need to resemble experimental recordings as closely as possible and contain the effects for which the analysis method searches. In our case this requires the simulation of a local cortical network representing roughly a cubic millimeter of cortex on which we can perform virtual experiments. This simulation is only briefly described here since this report is primarily about a new analysis tool for synfire chain activity.

The basic network architecture (Fig. 1) is the well-studied balanced random network model (Amit and Brunel 1997; Brunel 2000; Vreeswijk and Sompolinsky 1996). The dynamics of the network is adjusted to the asynchronous irregular (AI) regime to reproduce relevant aspects of cortical dynamics: low-rate firing of the individual neurons with Poisson process like interval statistics and large fluctuations of the membrane potential (see Fig. 1 for the full set of parameters). The model consists of 40,000 excitatory and 10,000 inhibitory integrate-and-fire point model neurons in the network.

The neurons establish a random number (drawn from a binomial distribution) of dendritic synapses with the excitatory population of neurons. The mean of the distribution is a fraction ε = 0.1 of the size of the excitatory population. This connection scheme permits multiple synapses between pairs of neurons. The scheme is also applied for establishing the dendritic synapses with the inhibitory population of neurons. The postsynaptic currents are described by alpha-functions with a rise time of 0.2 ms for the excitatory synapses and 0.6 ms for the inhibitory synapses. The excitatory postsynaptic potentials (EPSPs) have an amplitude of 0.1 mV; the inhibitory postsynaptic potentials (IPSPs) are (*g* = 6)-fold larger. The synaptic delays are randomly chosen from a uniform distribution between 0.5 and 3 ms. In addition to the explicitly modeled synaptic input all neurons receive a constant excitatory current of 350 pA describing the excitatory input from outside the local cortical volume. Thus the irregular activity is exclusively generated by the dynamics of the network and not imposed by an external source.

In contrast to the standard balanced random network architecture, in our network the excitatory-to-excitatory connections are excluded from the construction scheme described in the previous paragraph. This leaves the AI state largely unaffected because it is essentially generated by the interplay of the external excitatory suprathreshold drive and the local inhibition. The excitatory-to-excitatory subnetwork is, instead, purely composed of a superposition of *m* = 50 synfire chains. These chains are consecutively constructed. First, we randomly select *w* × *l* neurons (here: *w* = 100, *l* = 20) from the population of excitatory neurons without repetitions. In the next step, these neurons are connected into a feedforward subnetwork of *l* successive links of *w* neurons. Each neuron, except the neurons in the first link, establishes a dendritic synapse with all the neurons in the preceding group. This arrangement is called a complete divergent and convergent connectivity (Abeles 1982; see also Griffith and Horn 1963). The connection procedure is repeated *m* times. Thus an excitatory neuron may participate in several chains but occurs at most once in any given chain. The synaptic delays are drawn from the distribution specified earlier for the other synapses but is kept fixed for all synapses connecting two subsequent groups (links) of neurons. The amplitudes of postsynaptic potentials of synapses in the chains are fivefold larger than those of other excitatory connections.

In the absence of specific stimuli the synfire chains are not active and the excitatory neurons spike irregularly at a low rate. For the data presented herein we apply independent 1-Hz Poisson trains of stimuli to the synfire chains. A stimulus consists of a volley of 100 spikes. The individual spike times are drawn from a Gaussian distribution centered at the time of stimulation with SD of 1 ms. This volley is sent to all the neurons of the first group in a chain using synapses as strong as the intrachain synapses.

The dynamics are integrated on a discrete time grid of 0.1 ms using the exact integration method (Rotter and Diesmann 1999; see Morrison et al. 2007 for a detailed discussion of the application of the method to integrate-and-fire dynamics). We verified our results using a newly developed simulation technique capable of treating spike times and delays in continuous time in the framework of a globally time driven simulation scheme (Morrison et al. 2007). All simulations were carried out using the parallel computing capabilities (Plesser et al. 2007) of the NEST simulation tool (Gewaltig and Diesmann 2008) on an Infiniband compute cluster with 24 nodes with four compute cores and 8 GB of RAM per node.

The average spike rate in this network is 2.2 Hz for the excitatory neurons and 1 Hz for the inhibitory. The probability of synfire activity to survive until the end of the chain (>80% active neurons within 6 ms around the expected firing time) is 0.7–0.8.

Previous studies (Aviel et al. 2003; Mehring et al. 2003; see also the discussion in Vogels et al. 2005) have reported considerable difficulties in embedding even a single feedforward subnetwork capable of stable propagation of synchronous activity into the balanced random network. Our network differs from previous approaches in two essential properties. First, there is heterogeneity in the number of dendritic synapses. Second, the excitatory-to-excitatory subnetwork is purely composed of synfire chains. In fact, Tetzlaff et al. (2005) showed that the former condition is sufficient for stable propagation of a single synfire chain in an otherwise random balanced network.

However, in the present manuscript the network model is only employed as the generator of cortical spike data mimicking synfire activity. The details of the network dynamics and a parametric analysis of performance will be presented elsewhere.

### Data overview

Neuronal activity is commonly visualized as a raster display where each tick represents the spike activity of a particular neuron. In the neural net simulations we deliberately randomized the neuron number (identity) with respect to the embedded synfire structures to mimic the random sampling of real recordings. It is instructive, however, to remap some data so that adjacent lines in the raster are associated with particular links and chains. Thus (Fig. 2*A*) raster lines 8,001–8,100 represent the activity of neurons in the first link and lines 8,001–10,000 constitute the activity of the entire fifth (arbitrarily chosen) synfire chain. The strong, nearly vertical lines in the raster of Fig. 2*A* are firings of this synfire chain; other, background neuron firings in the raster represent either spontaneous activity or the activity of other chains, since each neuron in this chain is also a member of two to three other chains. Figure 2*B* shows a portion of this raster with high time resolution revealing details of several chain runs. The offset between short vertical segments corresponds to activations of successive links; the offset is the interlink propagation or synaptic delay. Note that both link synchrony and interlink propagation time are variable and that there is a failure of propagation in the middle run of the chain.

Without the remapping of neuron identities as in Fig. 2 synfire chain sequences are not visible in rasters. Figure 3*A* shows a 5-s segment of the activity of 150 excitatory neurons randomly selected (and not remapped according to chain membership) from a synfire simulation as described earlier (50,000 neurons, 50 fully connected synfire chains independently stimulated by low-frequency Poisson trains); the raster has no particular structure. The raster from a larger set of 10,000 excitatory neurons over the same time range is shown in Fig. 3*B*. Again, since the neurons are chosen randomly as they would be in a real recording, there are no structures like those in Fig. 2*B* that are directly ascribable to activations of individual synfire chains. However, there is an obvious broader vertical striping that connotes a coordinated fluctuation in the population activity; this striping becomes increasingly visible when the raster display shows a large number of neurons (contrast Fig. 3, *A* vs. *B*). It turns out that these coordinated rate fluctuations are the result of fluctuations in the number of synfire chains that have been activated at any time. Individual neurons, however, experience no rate increase. Figure 4 juxtaposes *1*) the times of chain stimulations with *2*) their (low-pass) temporal representation (i.e., the number of chains activated at any time) and *3*) the fluctuations of total population rate. These are clearly strongly, although not completely, related. We note that raster displays of simulation data from random nets, which contain no explicit synfire chains, also exhibit similar vertical striping (i.e., population rate fluctuations), but presumably caused by a different mechanism (Kriener et al. 2008).

### Analysis method

The basic data source is a large file specifying the firing times of each of many neurons. The exact format is a matter of laboratory preference and irrelevant here. For the work in this study networks were simulated with a temporal precision of 0.1 ms, the observed number of neurons *n* = 50,000, and the total observed time *T* = 100 s. The first step in the analysis is to define a bin size for time based on the expected interlink propagation delay (3 ms in Figs. 6–8 and 11–13; 2 ms in Fig. 10*C* for greater clarity). The number of bins is *K* = *T*/bin size, so that a bin is specified by the index *i* lying between 1 and *K*. For each such bin along the entire length of the data we find the set *S*(*i*) of all the neurons that fired in that (3-ms) interval. For large networks the *S*(*i*) can be of the order of ≥100. The size of each *S*(*i*) defines the population firing rate at the time of that bin. Figure 4 illustrates a smoothed version of *S*(*i*) size over a selected epoch.

The fundamental calculation we carry out is to compare the net activity at time bins *i* and *j.* This can be done by two conceptually different approaches. The first uses the *S*(*i*) sets and calculates the size of the intersection of *S*(*i*) and *S*( *j*)—i.e., the number of neuron identifications that appear in both sets. The appropriate normalization here is the size of the lesser of the two sets, thus producing a number between 0 and 1. With these definitions (1)

The second, more traditional approach, defines a 0/1 vector *V*(*i*) of length *N* at each bin *i*, where the entries of 1 indicate the identities of neurons active at that time. Now we take the dot product of the vectors and normalize by the product of their magnitudes (norms). This defines the cosine of the angle between the two vectors *V*(*i*) and *V*(*j*) and is a frequently used measure of similarity (or difference) between the vectors. Thus (2)

Note that the numerators of the two expressions are numerically precisely the same (for 0/1 vectors), but the normalizations make *M _{S}* generally larger and more nonlinear with respect to the number of neuron identity overlaps (higher contrast) than

*M*if the two sets

_{V}*S*(

*i*) and

*S*(

*j*) differ in size. Clearly

*M*(

*i*,

*i*) = 1 and

*M*(

*j*,

*i*) =

*M*(

*i*,

*j).*Therefore it suffices to inspect only the triangular matrix

*i*≤

*j*. With the stated bin structure,

*i*and

*j*are both discrete measures of time; we decided to locate the origin at the lower left corner for visualization as in the traditional two-dimensional correlation displays like the JPSTH (joint peristimulus time histogram; Aertsen et al. 1989). Here the diagonal of

*i*=

*j*extends from the lower left to the upper right and we are inspecting the lower triangular part of the intersection matrix. The programs have been designed to allow examination of any time ranges for

*i*and

*j*with

*i*<

*j*, allowing details in any portion of the triangular matrix to be visualized. For massively parallel recorded data with large

*N*the set approach is computationally less demanding and was used throughout this study. With either approach we plot

*M*(

*i*,

*j*) as a matrix, the color or gray level indicating the normalized activity intersection (or vector similarity) at the two time bins.

The measurement defined earlier for spike trains is a version of the “recurrence plot” of the mathematical literature (see Marwan et al. 2007 for a review). This was originally developed for following the states of a dynamical system and detecting approximately repeating states of the system. The formal definition of a recurrence plot involves a Heaviside function at a threshold value, so that matrix entries are limited to 0 or 1, and the emphasis is on detecting occurrence of states at times *i* and *j* that are more similar than the selected threshold value. In our application the state of the system is defined by the set of neurons active within a time window at time *i*, and we compute a matrix with continuous, not 0/1, values representing the degree to which the memberships of the active sets overlap at times *i* and *j.* Using adjustment of a gray or color scale we can emphasize either high or low values in the matrix.

### Implementation and run-time comparisons

All analysis and display programs were implemented in Matlab 2007a (The MathWorks) on a 3.1-GHz laboratory PC with 2-GB memory running the Windows XP system. Using this implementation we compared run times of the traditional pattern program (Abeles and Gerstein 1988) and the intersection matrix method for detecting and identifying members of synfire chains.

For run-time comparisons additional test data sets were made by choosing the activity of a random 500 or 100 neurons from the full 50,000 neuron network with 50 embedded synfire chains subsequently analyzed in Fig. 7. Data durations were 100 or 50 s. Except for the smaller number of neurons these test data had general characteristics similar to those described in *Data overview*.

The pattern program was used with a clock resolution of 2 ms and a time step of 2 ms. Repeating patterns <100 ms were accepted. Since many patterns arise by chance, a minimum of 21 iterations were required (original and 20 surrogates) to calculate significance. The intersection matrix program was used with a clock resolution of 0.1 ms, a 3-ms time bin, and window width of 1.5 s so that 66 or 33 time-stepped iterations were required to analyze the full data of 50 or 100 s. Direct significance testing is not necessary since, as subsequently demonstrated in *Control* and *Signal/noise, detection sensitivity*, the short diagonal lines in the matrix that are the signature of synfire chains do not arise by chance.

Timings are given in Table 1. Note that the time given for the pattern program does not include the as yet unknown time for sorting, collating, and comparing significant patterns to detect the meta structure and neuron identities (IDs) of synfire chains. In contrast the intersection matrix analysis is not only faster but directly detects synfire activity and requires only a simple additional computation to determine the IDs of the neurons in any chain. More detail is given in the next section.

### Identities of neurons in a synfire chain

The pattern programs sort spike patterns into groups according to complexity and repetition number and identify those groups whose count significantly exceeds corresponding counts found in surrogate data. All repeating patterns are consistent with an underlying synfire chain but could be generated by various other mechanisms. To resolve this issue it is necessary to examine various combinations of candidate patterns of various complexities and repetition number, searching for those that always interleave temporally in the same way (as would be expected from a synfire chain). Since there are very many candidate patterns (likely in the millions) this becomes a difficult combinatorial computation. We have not tried to write such programs as yet and can only guess at time estimates, but it is likely they would add considerably to the run time of the pattern programs.

By contrast, as shown in results, the intersection matrix programs directly identify the repeated runs of a synfire chain by short diagonal stripes, each of which represents a chain run repetition pair. A simple additional step is needed to actually identify the neurons that are involved in a given chain. Consider a chain that runs at times *t*_{1}, *t*_{2}, *t*_{3}, *t*_{4} … in the data set, with analysis as shown, for example, in Fig. 6. The short diagonal stripes representing the various run combinations are arranged in a vertical and horizontal array as in that figure. Contributions to each pixel located in any particular diagonal stripe represent the intersection of both chain and nonchain (accidental) neuron IDs active at those two times.

To identify the chain neurons we intersect pixel by pixel the identities in the diagonal stripe for *t*_{1}, *t*_{2} and the diagonal stripe for *t*_{1}, *t*_{3}, and then *t*_{2}, *t*_{3}, and so forth. The IDs of accidental firings at those times will disappear, in the multiple intersection process, because their activity is random with respect to the chain runs. The core IDs remaining in this process are precisely members of the chain. We have not yet written such programs, but their computational requirements are obviously small and would not impinge on the computation time comparisons of Table 1.

## RESULTS

### Synfire chain signature

##### TWO FIRINGS.

We first consider the matrix signature of a single synfire chain that fires twice, once beginning at (binned) time *t*_{1}, once at (binned) time *t*_{2}, as shown in the cartoons of Fig. 5. If for the moment we simplify by taking the interlink propagation time exactly equal to the chosen bin width, at *t*_{1} we observe a set *S*(*t*_{1}) of neurons that includes all active members of the first link in the chain. The set *S*(*t*_{2}) similarly includes all members of the first link in the chain that are active at *t*_{2}. Both sets can additionally include whatever other neurons in the net were active at their respective times. At *t*_{1} + 1 and *t*_{2} + 1 the two sets include all active members of the second link in each chain, and so on. Thus the matrix *M* defined earlier has high values (large intersection) at *t*_{1}, *t*_{2}; *t*_{1} + 1, *t*_{2} + 1; … *t*_{1} + *l*, *t*_{2} + *l*, where *l* is the chain length. A visualization of the matrix shows a short diagonal stripe of pixels with length *l* (∼20 in our cases) starting at time *t*_{1}, *t*_{2}. If for some reason activity does not run through the entire set of *l* links in one of the two chain firings, the diagonal stripe will be shorter as determined by the link at which propagation fails. If the link-to-link propagation time is not equal to the selected bin width, the short diagonal stripes will show some irregularity, which takes the form of shorter diagonal segments with 1-pixel offsets. (The situation is analogous to drawing an arbitrarily inclined line on a pixeled computer screen.) Finally, if the link-to-link propagation time is not constant—so that the two activations of the chain differ in detailed timing structure—there may also be gaps or deviations from 45° along the short diagonal stripe. All these distortions will depend in part on the relationship between selected bin width and the noisiness of the interlink delays. Thus a repeated synfire chain activation always produces a short approximately diagonal stripe in the intersection matrix, but its detailed appearance will be affected by minor moiré-like effects. Examples are shown in Figs. 6 and 8*A*.

##### MANY FIRINGS.

Now suppose our single synfire chain fires *k* > 2 times within our analysis time window. We now observe short diagonal stripes at *k* × (*k* − 1)/2 locations; if the successive firings are numbered, these will be at (binned) times 1, 2; 1, 3 … 1, *k*; 2, 3; 2, 4 … 2, *k*; … *k* − 2, *k* − 1; *k* − 2, *k*; *k* − 1, *k*. This situation is shown in Fig. 6 for *k* = 4. Finally, if there are several active chains, with each producing its appropriate set of short diagonal stripes, we obtain Fig. 7, which analyzes an arbitrarily chosen 1.5 s of the 100-s data set. Figure 7 is calculated from a full data set, in effect observing all its 50,000 neurons; the result clearly detects the existence and activity of some number of the embedded synfire chains. Although we have not done so in Fig. 7, we could obviously identify and label all diagonal stripes that originate from the activity of a particular synfire chain. Note that besides the prominent short diagonal stripes, Fig. 7 also shows some faint, nonuniform background that will be discussed later. Typical details of the short diagonal stripes together with chain labeling are shown in Figs. 6 and 8*A*, which are enlargements of portions of Fig. 7.

### Signal/noise, detection sensitivity

In the context of simulated data we may examine the performance limits of the *M _{S}* measure by reducing the number of observed neurons. The neurons are chosen randomly with no relation to the synfire chain structures. The analysis time window is constant for Fig. 8,

*A*–

*D*. When observing only 500 randomly chosen neurons out of the 40,000 excitatory network, the short diagonal stripes are somewhat noisier, so that not all elements of a particular stripe appear. However, comparing Fig. 8,

*A*and

*B*, most if not all stripes remain clearly visible. When only 200 randomly selected neurons are observed, we get Fig. 8

*C*. The stripes are now much more noisy, with many omissions, and not all stripes can be seen. This situation probably represents the current limit of sensitivity for detection of possible synfire chains in a data set. Such numbers, however, are just becoming possible in experiments, both those carried out by multiple electrodes or wires and with optical methods using Ca

^{2+}indicators.

When small numbers of neurons are observed, considerable visual enhancement for detecting the diagonal-stripe signature of synfire chain activity can be achieved by passing an appropriate short 45° diagonal filter (6–8 pixels along the diagonal, rectangular) over the noisy *M* matrix. When this is done for the data shown in Fig. 8*C* we get the matrix shown in Fig. 8*D*. All short diagonal stripes visible from the large data set (Fig. 8*A*) can now clearly be seen with the *N* = 200 analysis, so that detection of at least some synfire chains should be possible with even smaller *N*. Note that even with the noisy matrix of Fig. 8*C* the filtration has produced only a single spurious diagonal stripe (circled).

A more rigorous depiction of the signal/noise situation can be made by varying the number of observed neurons and comparing the average and SD of matrix values on and off a small diagonal stripe that had been detected with the full large data set. To have reasonable statistics we calculate over 40 disjoint sets of random choices of each number of neurons and in each iteration sum over 15 pixels along 45° both on and off the short diagonal stripe location. The results are shown in Fig. 9*A* and give a synfire chain detection threshold of about 50 observed neurons.

Another approach to the signal/noise separation does not depend on knowing about the location of short diagonal stripes, as mentioned earlier, and uses the entire matrix. Here we examine a histogram (probability density function [pdf]) of the pixel population as a function of filtered pixel value, comparing the results of both a 45 and a 135° narrow (6 pixel) filter. It is convenient to examine the results as a “survivor function,” i.e., (1 − cdf), where cdf is the cumulative distribution function. For the data shown in Fig. 8, *B* or *C* (*N* = 500 or 200) the survivor functions of the two filtered matrices all show a steep drop at a filtered pixel value 0.166 on a scale of 0–1 (Fig. 9*C*). This drop is an appropriate threshold value for comparing the pixel population in the two filter situations but is data and filter shape dependent, and must be examined for each application of such a calculation. Values larger than this threshold correspond to “signal” (the short diagonal stripes or equivalent pixel values in the 135° filtered case) and values below that threshold correspond to “noise” events. Again calculating over 40 disjoint sets of random choices of neuron at each number of neurons the comparison of pixel counts with values >0.166 after the 45 and 135° filters is shown in Fig. 9*B*. Thus our second type of signal/noise analysis yields a threshold for detecting synfire chain activity at about 100 observed neurons. Note that the validity of this second analysis does depend on using a matrix window that contains 5–10 of short diagonal stripes as in Fig. 8.

In this study we have used data from simulations of fully connected synfire chains; i.e., each neuron in a given link receives input from all neurons of the preceding link but with other parameter values set so there were about 30% failures of complete chain propagation. If we increase the link width but keep the number of converging inputs to each neuron the same as before, we have a diluted chain with partial connectivity. The detection of chain activity in such a case will be easier than that in the fully connected case. During a run of this effectively wider chain the set of active neurons in each link will be larger and the set of intersections between any two chain activations will also be larger (the numerator of *Eq. 1*). The denominator of *Eq. 1*, however, includes the amount of ongoing background activity at that time bin and thus increases by a smaller fraction than does the numerator. Therefore the normalized matrix pixel value corresponding to a particular link activity repetition will increase. One could also use just the unnormalized pixel value (i.e., just the numerator of *Eq. 1*), although variability in firing may make this measure more noisy. Other parametric chain changes such as reduction of the converging input number would cause many more propagation failures that in turn would shorten many of the signature diagonal lines in the matrix and might make detection of synfire activity more difficult.

In the signal-to-noise analyses of Fig. 9 we repeatedly sampled (without replacement) various numbers of excitatory neurons from a total pool of 50,000. Each excitatory neuron in the pool was involved in several of the 50 synfire chains and a given chain comprised about 4% of the total pool. Consider dilution of the data so that we leave intact the pool size and the neurons in some smaller number of chains, but replace the rest of the data by similar rate neurons that do not participate in chain activity. The only result in the intersection analysis matrix (like Fig. 7) will be a diminution of the number of short diagonal stripes, but not their intensity or contrast. If on the other hand the dilution process also replaces neurons in the remaining chains, the detection of synfire activity becomes more difficult and depends on the fraction of pool size represented by each remaining (now only partially observable) synfire chain. With a dilution such that the observable part of a given synfire chain constitutes about 1% of the neuron pool, detection requires simultaneous observation of 200–500 neurons (not shown).

### Intersecting synfire chains

In the case of intersecting synfire chains we might have some or all neurons in a given link of chain *A* also appear in a given link of chain *B*. In such a situation activity propagating along either chain could begin to propagate in the other chain, starting at the shared link. A cartoon of this situation is illustrated in Fig. 10*A*. Figure 10*B* is the associated raster with remapped identities showing activities of chain 1 in raster lines 1–2,000 and chain 2 in raster lines 2,001–4,000. Link 10 constitutes the same 100 neurons in both chains and it is clear that activity in chain 2 begins when link 10 of chain 1 is active (or vice versa if chain 2 is stimulated). In terms of the sets of active neurons in the appropriate time bins, this means after the 10th links are active that there will be an increase of active set size used as the denominator of *Eq. 1* (less than doubling since there is also background activity). There will also be a doubling of intersection value to a repetition of the (double) chain run (used as the numerator in *Eq. 1*). The increase in pixel value in the *M _{S}* matrix (

*Eq. 1*) is clearly visible in the gray levels of Fig. 10

*C*, which shows an enlarged and filtered image of the diagonal stripes corresponding to the three 2–1 chain-firing repetitions shown in Fig. 10

*B*. For this purpose we could also have used the unnormalized matrix of intersection values, i.e., the numerator of

*Eq. 1*alone; however, this may be more noisy if activities are variable between runs.

If instead link 10 of chain 1 were identical with link 5 of chain 2 we would expect a longer stripe with a profile along the stripe that increases and then returns to the single values. Thus within the limitations of count statistics it is possible in simple cases to infer details of intersecting synfire chains, using either normalized or unnormalized versions of the matrix.

More generally, detailed examination of the identities and numbers contributing to successive pixels of one or more short diagonal stripes could be used *1*) to determine the membership of a particular chain (described in more detail in the last section of methods) and *2*) to detect partial coupling between chains with resulting probabilistic ignition. Propagation failures would shorten the corresponding diagonal stripe and could be assigned to the responsible chain. However, general methods for appropriate comparisons of identity lists and for interpretation of such comparisons constitute a meta problem that remains to be addressed.

Another aspect of intersecting chains occurs in Fig. 8*A* in the short stripe labeled (1, 2). This stripe represents the neuron identity intersection of the second half of a full run of chain 2 with an earlier run of the second half of chain 2, which was caused by a full run of chain 1. Thus unlike most of the diagonal stripes, this does not represent a repetition of a directly stimulated synfire chain. Such a detailed explanation is possible here only because we know all about the data; it would not be possible in the case of experimental data.

### Background and blobs

The small diagonal lines in the intersection matrix that are the signature of synfire chain activity have pixel values near 0.6 and the gray scale of the figures in this study has been chosen to emphasize these larger values. As shown in Fig. 11, a change of gray scale in Fig. 7 allows examination of the many small values in the intersection matrix; these take the form of a general nonuniform background haze. There are broad, diffuse vertical and horizontal bands with blobs of somewhat higher density at their intersections. With high temporal resolution (not shown) individual pixels are spaced and have roughly the same low values, whether in or out of a broad band or blob; it is the density of pixels that creates these features. Referring back to Figs. 3*B* and 4, we note that the fluctuations of overall population firing rates have a similar time structure. During periods of higher population rate the set of active neurons in each time bin is larger and, correspondingly, the probability of a small identity overlap between two close times within a rate surge is higher, thus yielding a higher density of small pixels in the matrix. There is little overlap of neuron identities in the events leading to different broad bands or different blobs. We have demonstrated earlier (Fig. 4) that much of the population rate fluctuation in these data is caused by overlap of synfire chain stimulations. Similar population rate fluctuations in randomly connected nets without any synfire structures remain to be investigated by the active population intersection methods used herein.

### Controls

We use data from a number of simulated networks specified in detail in the following text. For each such data set the intersection matrix is featureless, with randomly distributed darker and lighter pixels. An example from a stimulated random network is shown in Fig. 12 . For all other control cases the intersection matrices are similar and are therefore not shown.

The first data set consists of 40,000 independent inhomogeneous spike trains, with a gamma interval distribution of order 4, and a noisy rate profile matching the average population firing rates in the synfire simulations. Individual neuron firing rates are partially correlated so that the population firing rate is a noisy oscillation similar to that shown in Fig. 4.

The second data set consists of spike trains from a randomly connected network, each neuron with the same general parameters and degree of inward and outward connectivity as in the synfire networks. Random activity is without specific external stimulation.

The third data set consists of spike trains from the same randomly connected network but with additional independent Poisson stimulation of 50 arbitrarily chosen groups of 100 neurons, matching the stimulation conditions used in the synfire simulations (Fig. 12).

The fourth data set consists of spike trains from a network with connectivity for 50 synfire chains but without the external stimulation (Y in Fig. 1) that was used to start chain activity. The network had only its intrinsic “spontaneous” activity.

In an additional control procedure we have increasingly disorganized the orderly timing progression of activity along a synfire chain. This was done by dithering—i.e., randomly shifting within some time window the roughly synchronous firing of all neurons in each link of each chain. For several values of dithering time window Fig. 13 shows a raster display of a typical run of dithered chain data and also shows a portion of the corresponding intersection matrix. Note that this procedure leaves intact the original approximate synchrony of activity in each 100-neuron link of the chain, but destroys the systematic sequential timing between successive links. Because the sequential nature of the link firing timing is increasingly disorganized (Fig. 13, *top* to *bottom*) the corresponding short diagonal lines in the matrix blur and dissolve. With a dither time window of 50 ms, which exceeds the approximately 45-ms duration of a normal synfire sequence, the intersection matrix has no remaining structure and shows only a random scattering of high-value pixels.

A final control calculation has to do with the possible consequences of bad spike sorting in multineuron recording. Spike-shape sorting is required when more than one neuron is recorded on a given (multi-) electrode. Although the literature is full of methods, sorting continues to be somewhat problematic and frequently results to some extent in multineuron rather than single-neuron spike trains (Harris et al. 2000). Spike train measures related to synchrony, auto- and cross-correlation, and time structure in general are considerably influenced by bad sorting (Gerstein 2000; Pazienti et al. 2006), although population measures are scarcely affected. It is therefore appropriate to examine the sensitivity of the intersection matrix to bad spike sorting. Data sets mimicking this situation were produced by overlaying sets of 5 or 10 randomly selected spike trains in the original 50,000 neuron synfire chain simulation file analyzed in Fig. 7, thus producing data sets with either 10,000 or 5,000 “mixed” spike trains. The resulting intersection matrices have somewhat larger values in the background gray pixels than in Fig. 7, but all the short diagonal lines continue to be dark and prominent. Thus even in this rather worst case of bad sorting scenario, the intersection matrix calculation is robust and continues to deliver the same information as in a clean, well-isolated neuron data situation.

These several control computations demonstrate clearly that the short diagonal structures seen in the intersection matrices of Figs. 6–8 and 11–13 do not arise accidentally from various other types of network activity and are, in fact, the unique signatures of the orderly and sequential activities of groups of neurons arranged in synfire chains.

## DISCUSSION

In this study we have presented a new method related to “recurrence plots” for detection and visualization of synfire chain activity. The calculation examines the normalized intersection of identities for neurons firing at any two distinct (binned) times *i* and *j*. With a binning resolution of milliseconds, the rapid and orderly sequence of link activation in a synfire chain produces a signature of short diagonal lines in the intersection matrix *M*(*i*, *j*).

We have shown that two related approaches are possible, one based on set intersections and one on vector representation of firing. Both types of measures are normalized to 0/1 but behave somewhat differently. The set intersection approach is computationally more convenient with very large numbers of neurons as in the simulated nets used herein. Compared with the vector approach, it also offers an improvement in contrast between high- and low-identity overlaps when the number of active neurons (set sizes) at times *i* and *j* differ. We have demonstrated that limited random samples of the order of 100 to 200 neurons from the large network are sufficient to detect synfire chains if they are present. This level of sensitivity should in the fairly near future allow experimental testing and either verification or rejection of the synfire idea. Current recording methods, electrode and optical, are allowing simultaneous observation of some 100 neurons in vivo, and rapid improvements can be expected.

We note that in the simulated data used here all excitatory neurons were involved in several synfire chains and that the neurons from any one synfire chain comprised about 4% of the total pool of observed neurons. Detection required a sample of about 100 neurons from the total pool. In experimental data from real nervous systems it may be that each synfire chain will constitute a smaller percentage of the observed neuron pool. This form of dilution would require larger random sampling numbers: if a single synfire chain constitutes only 1% of the total pool, it is necessary to sample 300–500 neurons for detection.

The vector representation of neuronal firing activity in some time window appears in the literature in a number of ways. It has been used in comparing simultaneously or sequentially recorded responses from several neurons to different and/or repeating stimuli (DiLorenzo 1989; Georgopoulos et al. 1983; Schild 1988). In these studies the response of each of the *N* neurons in the observed population is characterized by a spike count over some appropriate single time window for each repetition of the stimulus. Clustering in the *N*-dimensional population response space is used to predict the stimulus parameters from the neural activity. More detailed “activity vectors” (the time pattern of spikes after a stimulus) have been used to examine variability of single-unit responses to repeated stimulus presentations (Sasaki et al. 2007; Schreiber et al. 2003). Vector representation in the sense used here (an *N* length 0/1 vector representing data with millisecond time resolution) also underlies unitary-event analysis (Grün et al. 2002; Riehle et al. 1997). In that work such vectors are compared with what is expected from the firing rates of the individual neurons. Comparison of such vectors at two different times as in the present study (but with very low time resolution >1 s and, consequently, vectors that may not be 0/1) has been used to follow the slowly varying state of an observed neuronal network both by a matrix representation and with a principal components analysis (Sasaki et al. 2007). Although computationally similar in comparing vectors (or active neuron sets) at two times, in the present study we have used a time resolution of 3 ms, which is sufficiently small to make vectors that are always 0/1. Such time resolution is essential to detect the firing of synfire chains with the expected timescale of their ordered activity. A similar display has independently been developed in the field of music processing (Goto 2006).

With a time-bin resolution at the value of interlink delays expected in a synfire chain, the short diagonal stripes appearing in the intersection matrix each represent a particular repetition of activation of a particular chain. Repeated activity of a particular chain therefore appears as multiple short diagonal stripes arranged in rows and columns corresponding to the chain activation times. Different chains form separate series of short diagonal stripes.

Examination of the details in a particular short diagonal stripe yields considerable information about the synfire chain properties. For example, if in one repetition the activity of a synfire chain dies out before running through all links, all signature diagonal stripes involving that repetition will be shorter than the full chain length. If there are partially shared links between two chains this will be evident in the intersection identities and by modulation of the intersection numbers along the signature short diagonal stripe. Partial or full ignition of one chain by another can be extracted from the intersection identities. Link propagation times and their variation can be estimated from the way the signature diagonal is made up of short offset diagonal line segments.

The precision of synchronization in a synfire chain link (temporal jitter in spike coincidences) and the delay between the activation of successive neuron links are two distinct biophysically determined parameters in the simulation. Synchronization is mainly determined by the rise time of the postsynaptic potential (Gödeke and Diesmann 2008), whereas interlink delay is determined by conduction time or synaptic delay. Such parametric properties and issues will be examined in a subsequent study about the network simulation. In the present application, the bin size of our analysis needs to be adjusted to a compromise between the effects of these two parameters. When the synaptic delay is larger than the precision of synchronization and varies between synfire links, the signature of repeated synfire activity in the matrix remains a dark stripe at 45°. However, in this situation, the stripe may contain holes with values of background activity or may consist of several contiguous shorter diagonal stripes offset by one bin. For this reason processing the matrix with a 45° diagonal filter generally enhances the detection of synfire activity. Fluctuations in background activity may also scatter the matrix elements constituting the stripe away from the diagonal; a 45° diagonal filter extending a few bins perpendicular to the orientation of the diagonal may then be useful.

One may wonder why the simple analysis presented here may perform better than the elaborate statistical tools previously developed to detect the presence of individual spatiotemporal spike patterns. A conceptual difference is that our method automatically solves the puzzle of which patterns to consider as subpatterns of a larger one. We directly obtain a prominent diagonal matrix structure if synfire activity is present. Previously developed methods applied to such data detect many individual repeating patterns of neuron firings that do mostly originate from individual chains. However, such detected patterns must then be combined to determine the overall pattern of synfire activity, a computationally daunting project in a blind situation of real rather than simulated data. By intersecting identities of neurons firing at two different times our method integrates over the space of recorded neurons. By visual or filtered detection of diagonal structures that signify synfire activity, the method integrates over the temporal domain. This spatiotemporal integration provides sensitivity and robustness to the method.

## GRANTS

This work was partially supported by German-Israeli Project Cooperation Grant F1.2, European Union Fast Analog Computing with Emergent Transient States Grant 15879, and German Federal Ministry of Education and Research Grants 01GQ0413 and 01GQ0420 to Bernstein Center for Computational Neuroscience Berlin and Freiburg.

## Acknowledgments

We thank M. Abeles for urging development of a sensitivity test based on the full intersection matrix. Major parts of this work were conducted during a sequence of three scientific stays: August 2005, S. Grün and M. Diesmann with G. Gerstein at Jay, New York; September to December 2006, G. Gerstein with S. Schrader and M. Diesmann at Bernstein Center for Computational Neuroscience, Freiburg, Germany; and March 2007, G. Gerstein with S. Grün and M. Diesmann at RIKEN Brain Science Institute, Saitama, Japan.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society