With the ability to observe the activity from large numbers of neurons simultaneously using modern recording technologies, the chance to identify sub-networks involved in coordinated processing increases. Sequences of synchronous spike events (SSEs) constitute one type of such coordinated spiking that propagates activity in a temporally precise manner. The synfire chain was proposed as one potential model for such network processing. Previous work introduced a method for visualization of SSEs in massively parallel spike trains, based on an intersection matrix that contains in each entry the degree of overlap of active neurons in two corresponding time bins. Repeated SSEs are reflected in the matrix as diagonal structures of high overlap values. The method as such, however, leaves the task of identifying these diagonal structures to visual inspection rather than to a quantitative analysis. Here we present ASSET (Analysis of Sequences of Synchronous EvenTs), an improved, fully automated method which determines diagonal structures in the intersection matrix by a robust mathematical procedure. The method consists of a sequence of steps that i) assess which entries in the matrix potentially belong to a diagonal structure, ii) cluster these entries into individual diagonal structures and iii) determine the neurons composing the associated SSEs. We employ parallel point processes generated by stochastic simulations as test data to demonstrate the performance of the method under a wide range of realistic scenarios, including different types of non-stationarity of the spiking activity and different correlation structures. Finally, the ability of the method to discover SSEs is demonstrated on complex data from large network simulations with embedded synfire chains. Thus, ASSET represents an effective and efficient tool to analyze massively parallel spike data for temporal sequences of synchronous activity.
Synchronous input spikes to a receiving neuron are considered most effective in generating an output spike, as predicted by theoretical studies coining the term coincidence detector [1]. The argument rests on the premise that excitatory post-synaptic potentials (EPSPs) in the cortex are typically small in relation to the firing threshold, so that many EPSPs need to overlap to produce an output spike. Due to leak currents in the neuronal membrane, the firing threshold is reached with fewer spikes when these arrive synchronously at the post-synaptic neuron rather than sparsely, thus making the neuron behave like a coincidence detector. Experimental studies provide evidence for the existence of coincidence detectors (e. g. , [2]) and relate them to various mechanisms of spike-timing dependent plasticity [3,4] as well as to different encoding and decoding schemes [5,6]. Cortical anatomy supports such considerations. Individual neurons receive synaptic connections from a large number of neurons (on the order of 10,000 in the human cortex, see e. g. [7]) and project to a similar number of other cells. Such a connectivity structure combined with suitable synaptic delays may produce spatio-temporal spike patterns, as proposed in [8–11]. A simple case is represented by a temporal sequence of synchronous events (SSE), each event consisting of synchronous spikes from a group of neurons. The synfire chain model [8] is a neural network model that has been proposed to exhibit such activity through a suitable wiring of the neurons in the network in successive groups interconnected in a highly divergent and convergent manner. Each neuron of one group projects to several neurons of the next group, thus forming a chain structure. Under the assumption that the synaptic transmission delays from neurons of one group to neurons of the next group are identical, the synchronous stimulation of neurons in the first group leads to robust propagation of synchronous spiking activity through the chain [10] even in the presence of noise. The activation of a synfire chain would thus lead to the occurrence of an SSE. In order to assess whether SSEs are indeed observed in the brain and have a functional role, the spiking activity of several neurons needs to be recorded simultaneously, analyzed for temporal correlation and related to behavior. Under the assumption that an SSE occurs sparsely in a given data set, for instance in relation to a specific behavior, pairwise correlations between neurons engaged in the activity are weak and will not be detected by means of pairwise correlation analyses. For this reason, an analysis method was proposed in [12] that directly searches for SSEs in massively parallel spike data. The basic idea of the method presented in [12] is the following: after discretizing time in bins of a few ms (see Fig 1A), a synchronous event which repeats at two different time bins leads to a large number of neurons that are active in both bins. Building an intersection matrix I where each entry Iij represents the number of neurons active in both bins i and j, a group of synchronous events occurring in bins i and j results in a large overlap Iij compared to other entries of the matrix. An SSE which occurs twice produces in the matrix I a sequence of high-valued entries, which we name diagonal structures (DS), aligned parallel to the main diagonal (Fig 1B). A diagonal filter applied to the matrix I enhances the contrast between the DS and surrounding entries, mapping I into a filtered matrix F. The method was calibrated using data from large-scale synfire chain network simulations. When the neurons of a full chain or a large portion of it are observed, the method reveals the associated DS in the intersection matrix. However, when only a few hundreds of neurons (comparable to the number of cells that can be recorded simultaneously in vivo with modern electrophysiological techniques) are randomly sampled from the full network, the DSs are less visible and less likely to be continuous, making it impossible to isolate them from the surrounding entries. Indeed, this method leaves the judgment of whether individual entries are large or small, as well as the grouping of proximal entries into DSs, to visual inspection. As a consequence, the results are prone to subjectiveness and the procedure is not open to automation. In this paper we present a method, named ASSET (Analysis of Sequences of Synchronous EvenTs), which improves the approach proposed in [12] by providing a mathematical and fully automated detection of SSEs in parallel spike train data. The analysis features i) a statistical assessment of membership of individual entries of the intersection matrix to a DS, ii) a rigorous construction of individual DSs and associated SSEs by clustering, and iii) the reconstruction of the neuronal composition and occurrence times of each event in the found SSEs. The manuscript is organized as follows. In “Methods” we derive a statistical assessment of the membership of matrix entries to a DS based on two statistical tests for the significance of each individual entry and the joint significance of an entry and its neighbors, respectively. Entries passing both tests are then grouped together into individual DSs depending on their reciprocal distance by a clustering procedure. Once the DSs are identified, it is readily possible to reconstruct their neuronal composition. In “Results” we assess the performance and robustness of ASSET on various types of simulated data that replicate typical features of electrophysiological recordings. These include firing rate heterogeneity across neurons, variability over time and different types of correlation structure in the spiking activity. In addition, we demonstrate the ability of the ASSET analysis to find repeated SSE activity in simulated data generated by a neural network model of overlapping synfire chains as introduced in [12]. We conclude by discussing benefits and limitations of the ASSET method, open problems and future plans for the analysis of electrophysiological data.
Our approach builds on the notion of intersection matrix defined in [12]. Here we introduce this concept mathematically and provide definitions. Given a set of N parallel spike trains observed in the time interval [0, T], synchronous spike events across neurons are determined by discretizing the time interval into B adjacent time bins b1, b2, …, bB of identical width Δ = T/B (typically of a few ms), as illustrated in Fig 1A. Each set Si of spikes falling in time bin bi forms a synchronous event. The value |Si ∩ Sj| represents the number of neurons being active at both time bins bi and bj. The intersection matrix I is defined by I i j ≔ | S i ∩ S j | ∀ i, j = 1,2, …, B. In this setting I is a symmetric matrix, whose main diagonal contains the population histogram, i. e. the time histogram of the number of neurons simultaneously active in each time bin. A synchronous spike event occurring at two bins i and j, e. g. as a result of the repeated activation of the same synfire chain, results in a larger value for Iij compared to the chance level. A repeated SSE composed of lSSE successive synchronous events, each of which occurs at time bins (bir, bjr), r = 1,2…, lSSE, determines a sequence of large-valued entries Iir, jr in I—which we term a diagonal structure (DS) —as sketched in Fig 1B. To account for the variability of firing rates over time when comparing different entries of I, in [12] each entry of Iij was normalized by the number of neurons active in each bin bi and bj. After normalization, the entries take value 1 for a complete overlap and value 0 for no overlap. To enhance the contrast between high-valued entries belonging to the same diagonal structure and the surrounding entries, the matrix was further filtered by a linear kernel having an orientation parallel to the main diagonal. We take here a different approach, as outlined in the following section. To derive a measure of overlap that is independent of firing rates, we first need the probability mass function (pmf) pij (⋅) of each individual entry Iij in the intersection matrix, under the null hypothesis H0 that the spike trains under consideration are realizations of mutually independent Poisson processes. If the null hypothesis is rejected, i. e. if the observed value ξ taken by Iij is too large to be interpreted as chance, we classify the overlap Si ∩ Sj as a statistically significant repeated synchronous event. In “Results” we show that the statistics of the method are robust to deviations from Poissonianity as well as to the presence of various types of correlations other than repeated SSE activity. Under the stated hypotheses, the distribution of Iij is determined by the firing rate of each neuron k at the time bins bi and bj. Iij represents the (stochastic) number of neurons firing simultaneously in both bins and can thus be expressed as I i j = ∑ k = 1 N I i j (k), (1) where I i j (k) is a Bernoulli random variable taking value 1 if neuron k fires in both bins i and j, and 0 otherwise. The probability parameter p i j (k) of I i j (k) is related to the local firing rate λ i (k) of neuron k at bin bi and λ j (k) at bin bj by p i j (k) = (1 - e - λ i (k) Δ) (1 - e - λ j (k) Δ) for each i ≠ j, where each factor 1 - e - λ i (k) Δ is the probability for neuron k to emit at least one spike in the time bin i. The knowledge of the firing rate profiles of all neurons is therefore a prerequisite for the exact computation of the pmf pij (⋅) of Iij, i, j = 1, …, N. In light of Eq 1, Iij is a Poisson Bernoulli random variable [13]. Its pmf pij (⋅) is analytically given by p i j (ξ) = ∑ A ∈ P N; ξ ∏ k ∈ A p i j (k) ∏ h ∈ A C (1 - p i j (h) ), (2) where P N; ξ is the family of all possible subsets of ξ elements that can be extracted from the set {1, …, N}, A is one such subset and AC = {1, …, N}∖A is the complement of A. Thus pij (ξ) is a summation of (N ξ) addenda, for a total of 2N terms needed to compute pij (⋅). The computation is feasible for small N, but soon becomes prohibitive as N increases beyond a few dozens, as in the applications to large parallel recordings we are interested in. By use of Le Cam’s theorem [13], we approximate pij (⋅) by a Poisson density function p (λ) with rate parameter λ = ∑k λk p i j (ξ) ≃ p (λ) (ξ) = λ ξ e - λ ξ! . (3) The approximation error grows quadratically with the p i j (k) ’s and stays low if the p i j (k) ’s are sufficiently small, as shown in Fig 2. Given the firing rate profile (i. e. the rate at each time bin) of each neuron, we can thus calculate the pmf pij (⋅) and its cumulative distribution function (cdf) P i j (·) ≔ ∑ ξ < · p i j (ξ), either exactly from Eq 2 or approximately relying on Le Cam’s approximation in Eq 3. Transforming each entry Iij by its respective cdf P i j (·), we map the observed overlaps Iij to cumulative probabilities P i j ≔ P i j (I i j) = ∑ ξ < I i j λ ξ e - λ ξ! , (4) obtaining the probability matrix P ≔ (Pij) ij, as illustrated in Fig 3A. If the null hypothesis holds, then P i j (·) is the true probability distribution of the amount of overlap between bins bi and bj, and Iij is a realization from that probability distribution. If so, P i j: = P i j (I i j) takes N + 1 values x ∈ [0,1) (as many as the intersection values from 0 to N), and its cdf is the identity function over the domain of this set of values: Pr (Pij < x) = x for any x. If Pij is large (close to 1), the null hypothesis is rejected in favor of the alternative hypothesis that the observed overlap reflects active synchronization between the involved neurons at time bins bi and bj. After setting a significance threshold α1 (e. g. α1 = 0. 99), we classify all entries Iij for which Pij > α1 as statistically significant, along with the associated repeated synchronous events. Note that the sum in Eq 4 includes only values of ξ strictly lower than Iij. This choice, which assigns the observed value Iij to the critical region of the hypothesis test, ensures that 1 − Pij retains the property of a p-value, namely that Pr (1 − Pij ≤ y) = y for each y. 1 − Pij reflects the probability, computed under the null hypothesis, that Iij would take a value equal to or exceeding the observed value. We reject the null hypothesis if this probability is lower than 1 − α1 = 0. 01. A DS resulting from a repeated SSE, as sketched in Fig 1B, differs from the surrounding entries of the intersection matrix I due to not just one, but a sequence of entries with large values aligned around a diagonal. We devise a statistical test that exploits this feature to detect DSs in the intersection matrix. In the previous section we have already derived the probability matrix P as a transformation of I that normalizes raw intersection values by the local neuronal firing rates. We can now look for DSs in P rather than in I. Each entry Iij in the raw intersection matrix is tested for its individual significance and for the joint significance of its neighbors. If both tests pass, i. e. if Pij > α1 and Jij > α2, the entry is classified as belonging to one DS. Such entries are collected in a binary masked matrix M that takes value 1 if both tests pass, and 0 otherwise, M i j ≔ 1 { I i j > α 1 } · 1 { J i j > α 2 }, (8) as illustrated in Fig 3C. It remains to be established which entries belong together to the same DS. Intuitively, entries in the masked matrix that take value 1 belong to the same DS if close-by, and to different DSs if far apart. The masked matrix in Fig 3C shows for instance two clearly separated DSs. A suitable notion of “distance” for matrix entries should make entries falling in the same diagonal, i. e. aligned along the natural direction of a DS, closer together than entries aligned along the anti-diagonal. We introduce the following elliptic distance between any two matrix entries (i1, j1) and (i2, j2): d ρ ( (i 1, j 1), (i 2, j 2) ) ≔ 1 + (ρ - 1) · sin θ - π 4 · (i 2 - i 1) 2 + (j 2 - j 1) 2 / 2, where θ = arctan (j 2 - j 1 i 2 - i 1) is the angular coefficient of the line intersecting (i1, j1) and (i2, j2), the first square root factor is the Euclidean distance between the two points and ρ ≥ 1 is a stretching factor for angular coefficients deviating from π/4. d (⋅) grows as θ approaches 3π/4 or − π/4, i. e. the anti-diagonal orientation. For instance, dρ ( (i, j), (i+k, j+k) ) = k for any ρ, dρ ( (i, j), (i+k, j − k) ) = ρk and d ρ ( (i, j), (i, j + k) ) = [ 1 + 2 2 (ρ - 1) ] k. We set ρ = 5 and, based on the distance d5 (⋅, ⋅) of their positions in the matrix, group all entries Iij with Mij = 1 (see Eq 8) into clusters via a density based scanning (DBSCAN) algorithm [14]. The algorithm considers two entries as part of the same neighborhood if their distance is not larger than a maximum value ε. Neighborhoods sharing an entry are joined together and eventually classified as a cluster if they contain a minimum number l0 of entries. We set ε = 3. 5, thus allowing for a maximum of ⌊ε⌋ − 1 = 2 holes between two consecutive entries of a DS along the main diagonal, and l0 = 3, thus requiring a DS to reflect at least 3 repeated synchronous events. The elliptical neighborhood used when clustering should be contained into the kernel used to build the joint probability matrix J. The reason is that the first defines the “immediate” neighbors of the entry, while the second is meant to cover all entries that may belong to the same DS. Thus, the parameter ε should not be chosen larger than the kernel length lK, while the stretching coefficient ρ should be chosen such that the shorter axis of the ellipse fits into the kernel width wK. The values ε = 3. 5 and ρ = 5 we set for the validation satisfy these requirements. Fig 6 illustrates the matrix entries falling inside the ellipse (red dots) for various choices of the parameters ρ and ε, and the kernel (gray area) centered around the same entry. Entries in the matrix M not belonging to any cluster are discarded as events that do not reflect repeated SSE activity. Fig 3D shows the cluster matrix C assigning value 1 (colored in black) to entries belonging to a cluster, and 0 (white) to the others. The matrix contains two clusters composed of 5 entries each, corresponding to the 5 synchronous events highlighted in Fig 1A. Calculating the entries of the probability matrix as given by Eqs 2 and 3 requires the knowledge of the firing rate of each neuron over time. Firing rate estimation is a problem that has been targeted by a number of studies. The peri-stimulus time histogram (PSTH, [15]) is an estimate of the firing rate performed by discretizing time into adjacent bins and by counting the number of spikes falling into each bin. The bin width for the PSTH is typically larger (tens of ms) than the one used to define synchrony, as firing rates change on a slower time scale. The larger the bin width, the coarser but less biased the estimate. Normalizing the spike count in each bin by the bin width yields the rate of the process, i. e. the number of spikes per time unit. Kernel convolution [16] replaces each spike with a kernel (a probability density function) centered around the spike time and estimates the firing rate by the sum of these distributions. Formally, this is done by representing the spike time by a Dirac delta function centered around the spike time t* and by convolving it with the kernel. Following a standard choice, we specify the kernel as a normal distribution with assigned standard deviation σ, truncated at ±2. 7σ to yield a finite support. When setting the kernel width w* = 5. 4σ to a fixed value, we employ w* = 200 ms. Both PSTH and kernel convolution can be applied to the case when multiple independent, identically distributed trials of the activity of a neuron are available, by averaging the estimates obtained for each trial. For identical bin and kernel widths, the PSTH typically better represents sharp changes in the firing rate from one bin to the next, while kernel convolution yields smoother curves. Both estimates are parametric and require the choice of a bin- (kernel-) width. Methods have been recently proposed to determine the optimal bin- or kernel-width by minimizing the error between the true (unknown) rate and its estimate in some statistical sense (see [17,18]). These methods have been shown to outperform their fixed width variants and are particularly helpful when analyzing parallel spike train data from different neurons with different rates, where the optimal bin width varies across neurons. In the “Results” we compare the performance of ASSET employing either PSTH or kernel convolution or optimized-width kernel convolution [18] estimates of the firing rate profiles over an increasing number of trials.
To assess the quality of the method, we first measure its performance in the case when all assumptions entering the derivation of Eq 2 (alternatively Eq 3) and Eq 6 are met. In addition, we test the robustness of the statistics of the method with respect to deviations from these assumptions which are typically found in experimental data. In particular, we investigate how the following features of the data affect the performance of the method and test if these lead to false positive (FP) outcomes. The first four features relate to aspects of firing rates, such as various types of non-stationarities and rate correlations, i. e. correlations between the spike trains on a slower time scale than SSEs. The last two features relate to spike synchrony, however with a different organization than in SSEs. For the concrete test cases we formulate 10 different stochastic models of spiking activity (see Fig 8 as examples of the realizations), each including only one or a combination of the above-mentioned features, as summarized in Table 1. We use these models to generate background activity into which we subsequently embed the spiking activity corresponding to a repeated SSE. We provide here the definition of each model. The computational cost of the algorithm is almost entirely determined by the time required to evaluate the joint probabilities defined in Eq 6. The expression involves a nested sum of several terms. The number of terms grows with the number n of matrix entries covered by the kernel (determined by the kernel length lK and the kernel width wK) and with the number d of largest neighbors among which the joint significance is computed. These are free parameters of the analysis, whereas features of the data, such as the number of neurons or their firing rates, do not influence this step of the computation. As an estimate, for the values we employed in the manuscript (lK = 5, wK = 5, d = 5) the evaluation of a single entry took about 10 ms on a single core of a dual AMD 12-core Opteron 6174 machine with 64GB RAM using the Python code provided with this manuscript. Thus, the evaluation of a full matrix J of 200 × 200 entries took on average less than 7 minutes. However, the fact that single entries are evaluated independently may be easily exploited by parallelizing the analysis on multi-core machines or computer clusters, where each worker process is assigned to perform the computation for a subset of the matrix entries (see [37]).
Temporal sequences of synchronous spike events (SSEs) have been postulated as a working mechanism of activity propagation in the cortex [38–40]. The present manuscript introduces a novel statistical method for the detection of SSEs in massively parallel spike train data, named ASSET (Analysis of Synchronous Spike EvenTs). The method is inspired by a visual technique first proposed in [12], which represents the repeated occurrence of an SSE as a sequence of large entries along the diagonals of an intersection matrix (diagonal structure, or DS) that indicates for any two time bins the number of neurons firing in both bins. ASSET automatizes the detection process of the original visual technique, assesses the statistical significance of SSEs by exploiting the multiple evidence of its events to derive their joint significance, and determines the structure and neuronal composition of the identified SSEs. In evaluating the null distribution, the method accounts for the temporal profiles of the firing rates of the observed neurons. As such, it detects SSEs which cannot be explained on the basis of rate coding mechanisms, and thus arise from spike correlations on a shorter temporal scale. Rate correlation is understood as a conceptually different mechanism of information processing than spike synchrony (as for SSEs), because it corresponds to stochastic (probability <1) rather than reliable (probability 1) neuron activation [29,41–43]. We assessed the performance of ASSET in terms of false positive (FP) DSs found in various types of stochastic models which mimicked typical features of neuronal spike trains, such as variable firing rates in time or across neurons, different inter-spike interval distributions, and correlation structures differing from SSEs. We then additionally injected repeated SSE activity in the data to assess the power of the method in terms of true positive (TP) detections. The analysis performs two statistical tests on each pair of time bins, which amounts to tens of thousands of tests for a stochastic simulation of 1 s binned at 5 ms. In addition, entries passing the two tests have to lie close to each other in the matrix in order to be clustered into a common DS. To avoid FPs, the statistical threshold needs to be set to low values (here, 10−5). This is possible without incurring into large levels of FNs because the joint tail probabilities associated to the second test are very low (typically <10−5, and as low as 10−12) for the entries corresponding to the embedded SSEs. Indeed, the method shows high performance, i. e. FP and FN rates both close to 0. We did not need to further correct the statistical thresholds by the amount of tests performed (e. g. by Bonferroni or FDR correction), because already the set of combined requirements that need to be fulfilled to identify a DS (success of the two tests for matrix entries, and proximity of these entries in order to be clustered into a DS) makes the analysis very conservative. The underlying null hypothesis assumes independent Poisson spike trains, which enables an analytical formulation of the test statistics. The method proves to be robust to deviations from Poissonianity, such as a higher regularity of the inter-spike intervals, which may be observed in experimental data [22,44]. It is also selectively sensitive to SSEs, but not to other models of spike correlation, such as synchronous events not organized in a temporal sequence. Nevertheless, anticipating scenarios where strong correlations not forming an SSE might indeed bias the statistics, we proposed a Monte-Carlo approach to account for these correlations. To this end, we constructed the null hypothesis by estimating the probability to find a certain degree of pattern overlap from the repeated generation of surrogate intersection matrices. In our test data the Monte-Carlo approach yielded results comparable to the analytical approach, yet at a considerably higher computational cost. Furthermore, the method was able to distinguish SSEs from repeated precise temporal sequences of sharp, local rate transients from one group of neurons to another (rate propagation). Rate peaks increased the probability of the involved neurons to spike, which remained nevertheless a stochastic event, in contrast to SSE activity. Propagation of rate transients thus generated spike patterns that were different in composition (due to the stochastic activation of neurons) and less precisely timed than SSEs, but more closely resembled the latter as the rate modulation became higher and faster. It is likely that, for extremely high firing rates and very short rate jumps, ASSET would not distinguish the two models. The distinction between rate correlation and synchrony correlation (see e. g. [29]) has been formally questioned in [30]. Users of ASSET may want to identify waves of co-modulating rates like those defined in our model 6 as SSEs, rather than rejecting them. This is possible by decreasing the statistical thresholds α1 and α2 to less strict values than those used in this manuscript. ASSET critically relies on the estimation of the firing rate profile of each neuron to compute the expected overlap of neuron IDs in the intersection matrix and thus to estimate the significance of the observed overlap. Estimating firing rate profiles on the basis of single trial spike data typically requires some kind of sliding window approach, such as convolution of each spike with a temporal kernel [16]. Temporal averaging smears out peaks of the underlying original firing rate that occur on a shorter time scale than the window width, and creates artificial peaks if the window width is excessively short. Single-trial rate estimates obtained by kernel convolution in the presence of time-stationary firing rates yielded high performance of ASSET, but led to impaired performance when the rates were non-stationary in time on a fast time scale, in particular if the rate excursion was coherent across neurons. The reason is that smoothing by convolution underestimates positive rate peaks and thus the expected overlap, yielding FPs. Estimating the firing rates on the basis of trial averages solved the problem. We here tested three such approaches, namely the peri-stimulus time histogram (PSTH, [15]), a trial-averaged kernel convolution with fixed kernel width [16] and a trial-averaged kernel convolution with an optimized kernel width [18]. Already a small number of trials reduced the FP rate considerably for all three methods, although best performance was achieved for the optimized-width kernel convolution. Importantly, cross-trial rate estimation works under the assumption of identical rate profiles across trials. Deviations from this assumption lead to a wrong estimation of the rate in single trials, that is required to calculate the probability matrix, and thereby enhances the FPs. This bias is amplified if neurons exhibit cross-trial variability in a coherent manner [19]. Latency variability is a special instance of cross-trial non-stationarity which causes a mis-estimation around the rate onset. In some cases the onset variability of rates can be corrected for by choosing a more proper alignment of trials, e. g. to the stimulus or behavioral event related to the rate change [45]. If this is not possible, we suggest to generate the probability matrix P on the basis of surrogate spike data, e. g. by spike train shifting or spike dithering [20]. The details of such an approach, however, still need to be explored. We further investigated how the performance of ASSET relates to various other parameters of the SSEs, such as the number of its sequential synchronous events composing the SSE, the number of neurons in each synchronous event, and the total number of observed neurons. SSEs were statistically more significant and therefore easier to detect when they involved a larger fraction of the total neurons. However, they could be retrieved even when employing sub-optimal parameters such as a kernel length smaller than the length of the DS, given that the SSE involved enough neurons. In contrast, the method did not detect spatio-temporal patterns, i. e. a special case of SSE where each event is composed by a single spike only. Spatio-temporal patterns in a more general sense (with different time intervals between spikes) were suggested as signatures of synfire chain activity in data of low numbers of simultaneously recorded neurons [40,46,47]. A less constrained model of cortical processing (synfire braid) was proposed in [48], which incorporates synchronous input to individual neurons as in the synfire chain model, however transferred by connections of different temporal delays compensated by differences in their activation times. This model was further analyzed in [11] and termed polychronization. Since spike synchrony occurs with a temporal lag in this framework, one expects that, although ASSET is not designed to capture this type of coordinated spiking activity, it may still detect signatures of such activity by choosing a correspondingly larger bin width. We aim to explore such a scenario and other applications of ASSET to data from different types of network simulations that exhibit correlations on a fine temporal scale to study the potential of our method in identifying features of the underlying network model. Importantly, increasing the number of neurons to be analyzed does not increase the computational cost of the method. Indeed, in order to evaluate significance values, the method relies on expressions involving a mere sum of the firing rates of individual neurons, which are virtually instantaneous to evaluate. Rapid advances in electrophysiology will soon enable the simultaneous recording of thousands of neurons [49]. These data promise to expose concerted mechanisms of neuronal coding that remained invisible so far. Our method is designed to keep up with these advances, and to be applicable to the next generation of large-scale recordings of spike data. If an SSE occurs more than two times, it generates multiple DSs in the intersection matrix, each corresponding to one pair of occurrences. In [32] it was suggested to compute the overlap between triplets (or n-tuples) of bins rather than pairs, and generate a corresponding n-dimensional intersection matrix to visualize DSs in three (or n) dimensions. It is possible to extend this approach to ASSET and to exploit this higher-order evidence to increase the power of the method. This extension will be considered in future work. Finally, we demonstrated the efficacy of ASSET on data of large-scale simulations of a balanced random network with embedded synfire chains, which were previously generated and analyzed in [12] with the original visual method. ASSET fully reconstructed the synfire chains active in the considered time period and did not return additional FPs. Differently from the original technique in [12], SSEs were here determined on the basis of their statistical significance and the results were obtained by an automated analysis workflow. When analyzing real data, some parameters of the analysis such as the statistical thresholds α1 and α2 should be chosen optimally to minimize the risk of FPs, while at the same time not being excessively strict and thus missing true SSEs. Optimal values for the statistical thresholds can be inferred from independent surrogates of the original data, which can be created by one of several approaches, such as spike train shifting or spike dithering [50]. Such surrogate data contain slightly displaced spikes as compared to the original data, such that correlations in the original data (and in particular SSEs) are intentionally destroyed while other features of the data (e. g. firing rates or inter-spike interval regularities) are preserved. These uncorrelated data can be used to determine the expected value of entries in the probability and joint probability matrices under independence, and therefore to determine lower bounds for the thresholds α1 and α2 which ensure avoidance of FPs. Taking the least conservative of such values (the lower bounds) also minimizes the risk for FNs. Suitable values for other analysis parameters can be determined analogously. As illustrated in “Methods”, some of these parameters are tightly associated to putative features of the searched SSEs (for instance, the kernel length and kernel width to the temporal span and the wiggliness of the SSE, respectively), and may therefore be tuned in order to optimize the detection of SSE swith specific characteristics. Taken together, these results demonstrate that ASSET is a reliable and effective tool to detect and identify repeated sequences of synchronous spiking activity in massively parallel spike train data. Whether synchrony propagation constitutes a mechanism of information processing in the neuronal circuitry still remains an open question, that belongs to the more general debate about the role of fine temporal coding versus rate coding. Convincing theoretical arguments as well as experimental evidence have been provided for both processing mechanisms (see e. g. [2,6, 51] vs [26,52,53]). However, observing fine-scale temporal correlations requires the simultaneous observation and analysis of sufficiently large portions of the involved neuronal circuitry: the severe subsampling of such circuitry so far characterizing most available data prevents this analysis. Next-generation recording technology promises to expose this concerted activity, whose analysis may finally resolve the long-standing question about the role of millisecond-precise spike correlations in cognitive processing. Our work gives a contribution in this direction by providing, to the best of our knowledge, the first tool for a statistical analysis of synchrony propagation applicable to data of hundreds or thousands of simultaneously recorded neurons. The ASSET analysis method is available as part of the Electrophysiology Analysis Toolkit (Elephant; http: //neuralensemble. org/elephant/). In future work we plan to employ the method to investigate the presence of SSE activity in electrophysiological recordings from awake animals and study SSE occurrence in relation to behavior.
Neurons in the cerebral cortex are highly interconnected. However, the mechanisms of coordinated processing in the neuronal network are not yet understood. Theoretical studies have proposed synchronized electrical impulses (spikes) propagating between groups of nerve cells as a basis of cortical processing. Indeed, animal studies provide experimental evidence that spike synchronization occurs in relation to behavior. However, the observation of sequences of synchronous activity has not been reported so far, presumably due to two fundamental problems. First, the long-standing lack of simultaneous recordings of large populations of neurons, which only recently were enabled by advances in recording technology. Second, the absence of proper tools required to find these activity patterns in such high-dimensional data. Addressing the second issue, we introduce here a fully automatized mathematical method that advances an existing visual approach to identify sequences of synchronous events in large data sets. We demonstrate the efficacy of our method on a range of simulated test data that capture the characteristics and variability of experimental data. Our tool will serve future studies in their search for spike time coordination at millisecond precision in the brain.
action potentials
medicine and health sciences
neural networks
applied mathematics
membrane potential
electrophysiology
neuroscience
simulation and modeling
algorithms
probability distribution
mathematics
convolution
research and analysis methods
computer and information sciences
mathematical functions
animal cells
mathematical and statistical techniques
probability theory
cellular neuroscience
cell biology
kernel methods
neurons
physiology
biology and life sciences
cellular types
physical sciences
cognitive science
neurophysiology
ASSET: Analysis of Sequences of Synchronous Events in Massively Parallel Spike Trains