Leaders and followers: quantifying consistency in spatio-temporal propagation patterns

Most coincidence detectors rely on a coincidence window of fixed size τ [28, 29]. However, since in many cases it is very difficult to judge whether two spikes are coincident or not without taking the local context into account (see figure 2(a) for an example), [25] proposed a more flexible coincidence detection. This coincidence detection is scale- and thus parameter-free since the minimum time lag ${\tau }_{{ij}}^{(1,2)}$ at which two spikes ${t}_{i}^{(1)}$ and ${t}_{j}^{(2)}$ of spike trains $(1)$ and $(2)$ are no longer considered to be synchronous is adapted to the local firing rates according to

$$\begin{eqnarray}&&{\tau }_{{ij}}^{(1,2)}=\min \{{t}_{i+1}^{(1)}-{t}_{i}^{(1)},{t}_{i}^{(1)}-{t}_{i-1}^{(1)},{t}_{j+1}^{(2)}-{t}_{j}^{(2)},{t}_{j}^{(2)}-{t}_{j-1}^{(2)}\}/2.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

For some applications it might be appropriate to additionally introduce a maximum coincidence window ${\tau }_{\max }$ as a parameter. This way additional knowledge about the data (such as typical propagation speed) can be taken into account in order to guarantee that two coincident spikes are really part of the same propagation front. Such a maximum coincidence window will be used in the application to the El Niño climate data in section 3.3.

In normalization and windowing SPIKE-synchronization [27] has evolved so substantially from event synchronization that here we refrain from going into any detail on the original measure, but rather just mention the main improvements. For a thorough introduction to event synchronization please refer to the original paper [25], a more detailed comparison of the two measures can be found in [27].

The main difference is that SPIKE-synchronization [27] results in a discrete, not a continuous, spike-timing based profile. The coincidence criterion is quantified by means of a coincidence indicator

$$\begin{eqnarray}{C}_{i}^{(1,2)}=\left\{\begin{array}{ll}1 & \mathrm{if}\ {\min }_{j}(| {t}_{i}^{(1)}-{t}_{j}^{(2)}| )\lt {\tau }_{{ij}}^{(1,2)}\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

which assigns to each spike either a one or a zero depending on whether this spike is part of a coincidence or not. Note that here, unlike for event synchronization, the minimum function and the '<' guarantee that a spike can at most be coincident with one spike (the nearest one) in the other spike train. In case a spike is right in the middle between two spikes from the other spike train there is no ambiguity since this spike is not coincident with either one of them.

This unambiguity, illustrated in figure 2(b), is the essential property which allows the adaptive coincidence detection to act as a match-maker for the subsequent application of SPIKE-synchronization. Figure 2(c) shows examples, one with two coincident and one with two non-coincident spikes.

A multivariate version of SPIKE-synchronization can be defined by generalizing the bivariate coincidence detection of equation (2) to all pairs of spike trains (n, m) with $n,m=1,\ldots ,N$ and N denoting the number of spike trains:

$$\begin{eqnarray}{C}_{i}^{(n,m)}=\left\{\begin{array}{ll}1 & \mathrm{if}\ {\min }_{j}(| {t}_{i}^{(n)}-{t}_{j}^{(m)}| )\lt {\tau }_{{ij}}^{(n,m)}\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

Here ${\tau }_{{ij}}^{(n,m)}$ is defined equivalent to equation (1). Subsequently, for each spike of every spike train a normalized coincidence counter

$$\begin{eqnarray}&&{C}_{i}^{(n)}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{m\ne n}{C}_{i}^{(n,m)}\end{eqnarray} \tag{ 4 }$$

is obtained by averaging over all $N-1$ bivariate coincidence indicators involving the spike train n.

In order to obtain a single multivariate similarity profile we pool the spikes of all the spike trains as well as their coincidence counters:

$$\begin{eqnarray}&&\{{C}_{k}\}=\bigcup _{n}\{{C}_{i(k)}^{(n(k))}\},\end{eqnarray} \tag{ 5 }$$

where we map the spike train indices n and the spike indices i into a global spike index k denoted by the mapping i(k) and n(k).

Note that in case there exist perfectly coincident spikes, k counts over all of these spikes. From this discrete set of coincidence counters C_k the SPIKE-synchronization profile $C({t}_{k})$ is obtained via $C({t}_{k})={C}_{k}$. Finally, SPIKE-synchronization is defined as the average value of this profile

$$\begin{eqnarray}{S}_{C}=\left\{\begin{array}{ll}\displaystyle \frac{1}{M}{\displaystyle \sum }_{k=1}^{M}C({t}_{k}) & \mathrm{if}\ M\gt 0\\ \qquad 1 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

with $M={\sum }_{n=1}^{N}{M}^{(n)}$ denoting the total number of spikes in the pooled spike train.

This way we have used the same consistent framework for both the bivariate and the multivariate case. The former is just a special case of the latter. The interpretation is very intuitive: SPIKE-synchronization quantifies the overall fraction of coincidences. It reaches one if and only if each spike in every spike train has one matching spike in all the other spike trains (or if there are no spikes at all), and it attains the value zero if and only if the spike trains do not contain any coincidences. Examples for both of these extreme cases can be found in figures 3(a) and (b) and one intermediate example (random distribution of spikes among spike trains) is shown in figure 3(c). For a derivation of the expectation value for Poisson spike trains please refer to [30].

Zoom In Zoom Out Reset image size

In the multivariate analysis proposed in this paper, SPIKE-synchronization can be used to filter the input to the algorithm. In order to focus on propagation patterns within truly global events it is possible to set a threshold value C_thr for the SPIKE-synchronization profile $C({t}_{k})$. This way only spikes with a coincidence value higher than this parameter C_thr are taken into account, all the other noisy background spikes are simply ignored. This kind of filter will be used in the analysis of the neurophysiological datasets in section 3.2.

SPIKE-synchronization assigns to each spike of a given spike train pair a bivariate coincidence indicator. These coincidence indicators ${C}_{i}^{(n,m)}$, which are either 0 or 1, are then averaged over spike train pairs and converted into one overall profile $C({t}_{k})$ normalized between 0 and 1. In exactly the same manner SPIKE-order and spike train order assign bivariate order indicators to spikes. Also these two order indicators, the asymmetric ${D}_{i}^{(n,m)}$ and the symmetric ${E}_{i}^{(n,m)}$, which both can take the values −1, 0, or +1, are averaged over spike train pairs and converted into two overall profiles $D({t}_{k})$ and $E({t}_{k})$ which are normalized between −1 and 1. The SPIKE-order profile $D({t}_{k})$ distinguishes leading and following spikes, whereas the spike train order profile $E({t}_{k})$ provides information about the order of spike trains, i.e. it allows to sort spike trains from leaders to followers.

First of all, similar to the transition from the symmetric event synchronization to delay asymmetry, the symmetric coincidence indicator ${C}_{i}^{(n,m)}$ of SPIKE-synchronization (equation (3)) is replaced by the asymmetric SPIKE-order indicator

$$\begin{eqnarray}&&{D}_{i}^{(n,m)}={C}_{i}^{(n,m)}\cdot \mathrm{sign}({t}_{j^{\prime} }^{(m)}-{t}_{i}^{(n)}),\end{eqnarray} \tag{ 7 }$$

where the index $j^{\prime}$ is defined from the minimum in equation (2) as $j^{\prime} ={\arg \min }_{j}(| {t}_{i}^{(1)}-{t}_{j}^{(2)}| )$.

The corresponding value ${D}_{j^{\prime} }^{(m,n)}$ is obtained in an antisymmetric manner as

$$\begin{eqnarray}&&{D}_{j^{\prime} }^{(m,n)}={C}_{j^{\prime} }^{(m,n)}\cdot \mathrm{sign}({t}_{i}^{(n)}-{t}_{j^{\prime} }^{(m)})=-{D}_{i}^{(n,m)}.\end{eqnarray} \tag{ 8 }$$

Therefore, this indicator assigns to each spike either a 1 or a −1 depending on whether the respective spike is leading or following a coincident spike from the other spike train. The value 0 is obtained for cases in which there is no coincident spike in the other spike train (${C}_{i}^{(n,m)}=0$), but also in cases in which the times of the two coincident spikes are absolutely identical (${t}_{j^{\prime} }^{(m)}={t}_{i}^{(n)}$).

The multivariate profile $D({t}_{k})$ obtained analogously to equation (5) is normalized between 1 and −1 and the extreme values are obtained if a spike is either leading (+1) or following (−1) coincident spikes in all other spike trains. It can be 0 either if a spike is not part of any coincidences or if it leads exactly as many spikes from other spike trains in coincidences as it follows. From the definition in equations (7) and (8) it follows immediately that C_k is an upper bound for the absolute value $| {D}_{k}|$.

While the SPIKE-order profile can be very useful for color-coding and visualizing local spike leaders and followers (figure 4(a)), it is not useful as an overall indicator of spike train order (figure 4(b)). The profile is invariant under exchange of spike trains, i.e. it looks the same for all events no matter what the order of the firing is (in our example only the last event looks slightly different since one spike is missing). Moreover, summing over all profile values, which is equivalent to summing over all coincidences, necessarily leads to an average value of 0, since for every leading spike (+1) there has to be a following spike (−1).

Zoom In Zoom Out Reset image size

So in order to quantify any kind of leader-follower information between spike trains we need a second kind of order indicator. The spike train order indicator is similar to the SPIKE-order indicator defined in equations (7) and (8) but with two important differences. Both spikes are assigned the same value and this value now depends on the order of the spike trains:

$$\begin{eqnarray}&&{E}_{i}^{(n,m)}={C}_{i}^{(n,m)}\cdot \left\{\begin{array}{l}\mathrm{sign}({t}_{j^{\prime} }^{(m)}-{t}_{i}^{(n)})\quad \mathrm{if}\quad n\lt m\\ \mathrm{sign}({t}_{i}^{(n)}-{t}_{j^{\prime} }^{(m)})\quad \mathrm{if}\quad n\gt m\end{array}\right.\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{E}_{j^{\prime} }^{(m,n)}={E}_{i}^{(n,m)}.\end{eqnarray} \tag{ 10 }$$

This symmetric indicator assigns to both spikes a +1 in case the two spikes are in the correct order, i.e. the spike from the spike train with the lower spike train index is leading the coincidence, and a −1 in the opposite case. Once more the value 0 is obtained when there is no coincident spike in the other spike train or when the two coincident spikes are absolutely identical.

The multivariate profile $E({t}_{k})$, again obtained similarly to equation (5), is also normalized between 1 and −1 and the extreme values are obtained for a coincident event covering all spike trains with all spikes emitted in the order from first (last) to last (first) spike train, respectively (see the first two and the last four events in figure 4). It can be 0 either if a spike is not a part of any coincidences or if the order is such that correctly and incorrectly ordered spike train pairs cancel each other. Again, C_k is an upper bound for the absolute value of E_k.

In contrast to the SPIKE-order profile D_k, for the spike train order profile E_k it does make sense to define an average value, which we term the synfire indicator:

$$\begin{eqnarray}&&F=\displaystyle \frac{1}{M}\displaystyle \sum _{k=1}^{M}E({t}_{k}).\end{eqnarray} \tag{ 11 }$$

The interpretation is very intuitive. The synfire indicator F quantifies to what degree the spike trains in their current order resemble a perfect synfire pattern. It is normalized between 1 and −1 and attains the value 1 (−1) if the spike trains in their current order form a perfect (inverse) synfire pattern. This means that all spikes are coincident with spikes in all other spike trains and that all orders from leading (following) to following (leading) spike consistently reflect the order of the spike trains. The synfire indicator is 0 either if the spike trains do not contain any coincidences at all or if among all spike trains there is a complete symmetry between leading and following spikes.

The spike train order profile $E({t}_{k})$ for our example is shown in figure 4(c). In this case the order of spikes within an event clearly matters. The synfire indicator F is slightly negative indicating that the current order of the spike trains is actually closer to an inverse synfire pattern.

Given a set of spike trains we now would like to sort the spike trains from leader to follower such that the set comes as close as possible to a synfire pattern. To do so we have to maximize the overall number of correctly ordered coincidences and this is equivalent to maximizing the synfire indicator F. However, it would be very difficult to achieve this maximization by means of the multivariate profile $E({t}_{k})$. Clearly, it is more efficient to sort the spike trains based on a pairwise analysis of the spike trains. The most intuitive way is to use the anti-symmetric cumulative SPIKE-order matrix

$$\begin{eqnarray}&&{D}^{(n,m)}=\displaystyle \sum _{i}{D}_{i}^{(n,m)}\end{eqnarray} \tag{ 12 }$$

which sums up orders of coincidences from the respective pair of spike trains only and quantifies how much spike train n is leading spike train m (figure 5).

Zoom In Zoom Out Reset image size

Hence if ${D}^{(n,m)}\gt 0$ spike train n is leading m, while ${D}^{(n,m)}\lt 0$ means m is leading n. If the current spike train order is consistent with the synfire property, we thus expect that ${D}^{(n,m)}\gt 0$ for $n\lt m$ and ${D}^{(n,m)}\lt 0$ for $n\gt m$. Therefore, we construct the overall SPIKE-order as

$$\begin{eqnarray}&&{D}_{\lt }=\displaystyle \sum _{n\lt m}{D}^{(n,m)},\end{eqnarray} \tag{ 13 }$$

i.e. the sum over the upper right tridiagonal part of the matrix ${D}^{(n,m)}$.

After normalizing by the overall number of possible coincidences, we arrive at a second more practical definition of the synfire indicator:

$$\begin{eqnarray}&&F=\displaystyle \frac{2{D}_{\lt }}{(N-1)M}.\end{eqnarray} \tag{ 14 }$$

The value is identical to the one of equation (11), only the temporal and the spatial summation of coincidences (i.e., over the profile and over spike train pairs) are performed in the opposite order.

Having such a quantification depending on the order of spike trains, we can introduce a new ordering in terms of the spike train index permutation $\varphi (n)$. The overall synfire indicator for this permutation is then denoted as F_φ. Accordingly, for the initial (unsorted) order of spike trains ${\varphi }_{u}$ the synfire indicator is denoted as ${F}_{u}={F}_{{\varphi }_{u}}$.

The aim of the analysis is now to find the optimal (sorted) order ${\varphi }_{s}$ as the one resulting in the maximal overall synfire indicator ${F}_{s}={F}_{{\varphi }_{s}}$:

$$\begin{eqnarray}&&{\varphi }_{s}:{F}_{{\varphi }_{s}}=\mathop{\max }\limits_{\varphi }\{{F}_{\varphi }\}={F}_{s}.\end{eqnarray} \tag{ 15 }$$

This synfire indicator for the sorted spike trains quantifies how close spike trains can be sorted to resemble a synfire pattern, i.e., to what extent coinciding spike pairs with correct order prevail over coinciding spike pairs with incorrect order. Unlike the synfire indicator for the unsorted spike trains F_u, the optimized synfire indicator F_s can only attain values between 0 and 1 (any order that yields a negative result could simply be reversed in order to obtain the same positive value). For a perfect synfire pattern we obtain F_s = 1, while sufficiently long Poisson spike trains without any synfire structure yield ${F}_{s}\approx 0$.

The complexity of the problem to find the optimal spike train order is similar to the well-known traveling salesman problem [31]. For N spike trains there are $N!$ permutations φ, so for large numbers of spike trains finding the optimal spike train order ${\varphi }_{s}$ is a non-trivial problem and brute-force methods such as calculating the F_φ-value for all possible permutations are not feasible. Instead, one has to make use of methods such as parallel tempering [32] or simulated annealing [33] to search for the optimal order. Here we choose simulated annealing, a probabilistic technique which approximates the global optimum of a given function in a large search space. In our case this function is the synfire indicator F_φ (which we would like to maximize) and the search space is the permutation space of all spike trains. We start with the F_u-value from the unsorted permutation and then visit nearby permutations using the fundamental move of exchanging two neighboring spike trains within the current permutation. The update of the synfire indicator when exchanging the spike trains k and $k+1$ is simply given by ${\rm{\Delta }}F=-2{D}^{(k,k+1)}$. All moves with positive ${\rm{\Delta }}F$ are accepted while the likelihood of accepting moves with negative ${\rm{\Delta }}F$ is decreased along the way according to a standard slow cooling scheme. The procedure is repeated iteratively until the order of the spike trains no longer changes or until a predefined end temperature is reached.

In figure 6 we show the complete SPIKE-order analysis including the results for the sorted spike trains. The sorting of the spike trains maximizes the synfire indicator as reflected by both the normalized sum of the upper right half of the pairwise cumulative SPIKE-order matrix (equation (14), figure 6(c)) and the average value of the spike train order profile $E({t}_{k})$ (equation (11), figure 6(d)). Finally, the sorted spike trains in figure 6(e) are now ordered such that the first spike trains have predominantly high values (red) and the last spike trains predominantly low values (blue) of $D({t}_{k})$.

Zoom In Zoom Out Reset image size

The complete analysis returns results consisting of several levels of information. Time-resolved (local) information is represented in the spike-coloring and in the profiles D and E. The pairwise information in the SPIKE-order matrix reflects the leader-follower relationship between two spike trains at a time. The synfire indicator F characterizes the closeness of the dataset as a whole to a synfire pattern, both for the unsorted (F_u) and for the sorted (F_s) spike trains. Finally, the sorted order of the spike trains is a very important result in itself since it identifies the leading and the following spike trains.

As a last step in the analysis we evaluate the statistical significance of the optimized synfire indicator F_s. What we would like to estimate is the likelihood that for the given total number of coincidences the prevalence of correctly ordered spike pairs (as quantified by the optimized synfire indicator) could have been obtained by chance. If all coincident spike pairs would be independent, the probability distribution would be strictly binomial and we could calculate this likelihood analytically. However, the pairwise spike orders in coincident events involving multiple spike trains are not independent from each other, and so instead we estimate the likelihood numerically using a set of carefully constructed spike order surrogates.

For each surrogate (figure 7(a)) we maintain the coincidence structure of the original spike trains by preserving the SPIKE-synchronization values of every individual spike. However, we destroy the spike order patterns by swapping the order of the two spikes in a sufficient number of randomly selected coincident spike pairs. Note that the generation of surrogates takes place not on the level of spike times but on the level of order values (the x-axis in figure 7(a) is labeled 'time index', not 'time'). Spike trains with swapped spike times would have different interspike intervals, and this would alter the results of the coincidence criterion in equation (1) and change the value of SPIKE-synchronization. This in turn would make the desired evaluation of pure spike order effects difficult.

Zoom In Zoom Out Reset image size

In the implementation, from one spike order surrogate to the next the number of spike order swaps is set to the number of coincident spikes in the spike train set, such that all possible spike order patterns can be reached. Only for the first surrogate, since it starts from the original spike trains, we swap twice as many coincidences in order to account for transients. After each swap we take extra care that all other spike orders that are affected by the swap are updated as well. For example, if a swap changes the order between the first and the third spike in an ordered sequence of three spikes, we also swap both the order between the first and the second as well as the order between the second and the third spike.

For each spike train surrogate we repeat exactly the same optimization procedure in the spike train permutation space that is done for the original dataset. The original synfire indicator is deemed significant if it is higher than the synfire indicator obtained for all of the surrogate datasets (this case will be marked by two asterisks). Here we use s = 19 surrogates for a significance level of ${p}^{* }=1/(s+1)=0.05$. As a second indicator we state the z-score, e.g., the deviation of the original value x from the mean μ of the surrogates in units of their standard deviation σ:

$$\begin{eqnarray}&&z=\displaystyle \frac{x-\mu }{\sigma }.\end{eqnarray} \tag{ 16 }$$

Results of the significance analysis for our standard example are shown in the histogram in figure 7(b). In this case the absolute value of the z-score is smaller than one and the p-value is larger than p^* and the result is thus judged as statistically non-significant.

In case the initial sorting of the spike trains is used to test a specific hypothesis there also exists a straightforward procedure to test the statistical significance of the synfire indicator F_u for the unsorted spike trains. In this case no optimization of the synfire indicator is required, rather the synfire indicator F_u for the initial sorting is compared against synfire indicators obtained for random permutations of the spike trains. This kind of significance test will be used in section 3.2.

We start with examples covering the two extreme cases of a perfect synfire pattern and a completely random spike train set. First, in figure 8 we apply the algorithm to a perfect inverse synfire pattern for which the spike trains are initially sorted from follower to leader. Therefore, the synfire indicator of the unsorted spike trains yields its minimum value of ${F}_{u}=-1$. Sorting just reverses the order of the spike trains and in consequence the maximum value of F_s = 1 is obtained. Any shuffling of spike orders necessarily destroys the synfire pattern and thus leads to much lower values of the synfire indicator. Accordingly, the surrogate test (figure 8(f)) shows that the statistical significance of the original synfire indicator is very high.

Zoom In Zoom Out Reset image size

The other extreme case is Poisson spike trains (figure 9) for which the arrival times of spikes are completely random and without any preferred order. For this realization the synfire indicator F_u for the unsorted spike trains happens to be slightly negative indicating that the spike trains are closer to an inverse synfire pattern than to a synfire pattern. The absolute value F_s after sorting is higher. The fact that both of these values are non-zero is due to the finite size effect caused by the limited number of spikes. For more and more spike trains and/or spikes the expectation value even for the sorted case would converge towards zero. As expected, the surrogate test shows that the order for the original spike trains is not statistically distinct from the order of the surrogate spike trains (there is no preferred order that can be destroyed by the shuffling) and, accordingly, the value of the original synfire indicator is revealed to be clearly non-significant.

Zoom In Zoom Out Reset image size

The third example in figure 10 shows a mixture of these two extremes, Poisson spike train interspersed with spike trains that contain a perfect inverse synfire pattern (plus random spikes). Sorting the spike trains restores the correct order of the synfire pattern spike trains within the Poisson spike train. The synfire indicator for the sorted spike trains F_s for this mixed example is actually almost identical to the value obtained for the Poisson spike trains in figure 9, but this time the surrogate test reveals the value to be highly significant. These two examples combined illustrate nicely that the synfire indicator and the surrogate analysis provide complementary information. In the mixture example of figure 10 there are many more random Poisson spikes than ordered synfire pattern spikes. According to the synfire indicator, these two types of spikes together appear to be as ordered as the spikes of the shorter but purely random Poisson spike trains in figure 9. However, the synfire indicator is strongly influenced by the statistics of the dataset and thus is in itself not sufficient to reliably compare two datasets with widely different number of spike trains and spikes. The surrogate analysis, on the other hand, can be used to compare datasets of different size since by preserving the spike numbers in the surrogates it explicitly takes the statistics of each dataset into account.

Zoom In Zoom Out Reset image size

In order to apply the spike train order algorithm to real neurophysiological data, we analyzed data recorded via fast multicellular calcium imaging in acute CA3 hippocampal brain slices from juvenile mice. In the juvenile hippocampus, the CA3 region is the origin of a stereotypical network phenomenon of wavelike propagating activity termed giant depolarizing potentials (GDPs [34]). In previous studies, GDPs have been used to investigate the topology of networks and the role of hub cells [23] as well as to reveal the deterministic and stochastic processes underlying spontaneous, synchronous network bursts [35]. Due to the distinct architecture and the repetitive nature of the GDPs this experimental setup offers a very suitable test case for our synfire pattern analysis (for more background and a detailed description of the experimental methods refer to appendix A).

The first dataset analyzed in figure 11 includes 13 GDPs over a bit more than 6 min. Almost all GDPs involve the whole network. Here, as for all other neurophysiological datasets analyzed, initially the spike trains are sorted according to their firing rate such that the sparsely spiking neurons are on top and the most active neurons at the bottom (figure 11(a)). This specific sorting allows us to test the hypothesis that the neurons which fire almost exclusively within the GDPs and are very sparse on background activity might have a stronger role in initiating GDPs and tend to lead, whereas the more regularly spiking neurons might tend to follow. If this would be the case one would expect a very high value for the initial synfire indicator F_u. However, according to figure 11(b) the actual value is very close to zero and actually slightly negative. A statistical significance test using random permutations of spike trains (see section 2.5) indeed proves the synfire indicator of the unsorted spike trains F_u to be non-significant (result not shown). A further indicator for this is the fact that the order of the sorted spike trains is very different from the initial order, as can be seen by comparing the color bars on the left of figures 11(a) and (e). The color-coding of the GDPs exhibits typically a slightly noisy transition from leader (red) to follower (blue). The synfire indicator for the sorted spike trains F_s is also much higher (figure 11(d)). Finally, the surrogate analysis (figure 11(f)) shows this result to be highly significant.

Zoom In Zoom Out Reset image size

However, the spiking in figure 11 consists not only of the GDPs. Most neurons exhibit at least to some extent spontaneous background activity, the ones at the top of the initial sorting less than the ones at the bottom. The spikes in this background activity are typically coincident with only few other spikes and do not take part in any propagation patterns (note their green color which indicates SPIKE-order values close to zero). So in the context of our synfire pattern analysis this is just noise that leads to a decrease of the synfire indicator. There is a straightforward way to disregard these background spikes by setting a threshold value C_thr for the SPIKE-synchronization profile $C({t}_{k})$. Only spikes with a coincidence value higher than C_thr are taken into account, all other spikes are simply ignored. The result of this background correction can be seen in figure 12 for the same dataset already used in figure 11. Focusing the analysis on the reliable GDPs leads to an increase of the synfire indicator from 0.284 to 0.438.

Zoom In Zoom Out Reset image size

As already mentioned before, one of the main results of our analysis is the sorted order of the spike trains itself. For these neurophysiological data it allows to identify the leading and the following neurons in the network and to project this information back on the recording setup. This is shown in figure 13 where we have color-coded the optimized spike train order obtained in figure 12 within a 2D-plot of the neurons recorded from the hippocampal slice. For this example there appears to be a clear overall propagation from right to left but there is also a considerable degree of variability which might be due to a non-trivial connectivity within the network.

Zoom In Zoom Out Reset image size

In figure 14 we apply the SPIKE-order analysis to a second dataset recorded from a different slice, again focusing on the order within the global events only. Here we also added one new feature, the mean value of the spike train order $E({t}_{k})$ for each global event (we use the maxima and minima of the SPIKE order profile $D({t}_{k})$ to delineate the GDPs). This again emphasizes the time-resolved nature of the SPIKE order and the spike train order indicators.

Zoom In Zoom Out Reset image size

Overall, we have analyzed neurophysiological datasets from four hippocampal slices exhibiting an average of 7.75 GDPs. We obtained an average value for SPIKE-synchronization of 0.59 before focusing on the GDPs (as in figure 11) and 0.92 after (as in figure 12). With or without this focus the synfire indicator for the initial spike train sorting F_u was very close to zero and in all cases proved to be non-significant when tested against synfire indicators obtained for random permutations of the spike train order. Since the initial sorting was based on overall firing rate of the neurons, this signifies that the hypothesis that the low-firing neurons which are basically only active during the GDPs might have a stronger role in initiating GDPs can be rejected. For the sorted spike trains the synfire indicator F_s was 0.20 for all spikes and 0.42 for the spikes within the GDPs only. Suppressing the effect of the noisy background spikes in the analysis thus leads to an average increase of the synfire indicator by about $110 \%$. Finally, according to the surrogate analysis described in section 2.5 the synfire indicator for the sorted spike trains F_s yielded a statistically significant result for all datasets analyzed.

So overall we can conclude that the GDPs recorded in brain slices from juvenile mice are distinguished by a very high consistency of their spatio-temporal propagation patterns. However, it is interesting to note that this consistency does not hold when comparing different slices. In the datasets analyzed in this paper we find examples of both propagation in the direction of CA2 as well as propagation towards the dentate gyrus. This is consistent with results reported in [35].

Although being developed mainly for neuroscientific data (spike train recordings), the spike train order approach presented in this paper can be applied in many other contexts as well. One particular field where event-based analysis is employed very frequently is climate science, see e.g. [20, 36].

In the following we will use the new spike train order method to analyze the sea surface temperature (SST) in the central and eastern tropical Pacific Ocean to identify the propagation patterns connected with the warm phase of the El Niño Southern Oscillation (ENSO). Predicting the occurrence and strength of El Niños is very important due to the vast ecological and economical effects [37] and therefore this phenomenon has been studied extensively in terms of both intensity (e.g. [38, 39]) and frequency (e.g. [40, 41]). However, here we will not discuss El Niño prediction, but focus on the longitudinal propagation pattern of the El Niño events within the Pacific Ocean [24, 42]. The region used here to analyze the SST is depicted in figure 15 (dashed line), and an examplary smoothed time profile for the center of the analyzed region (215 °E) is shown in figure 16. Further details on the region, the dataset and data preprocessing (smoothing) is outlined in appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Starting from the smoothed SST profiles, we performed an event detection by identifying the onset of the SST anomaly connected with an El Niño in terms of an upward crossing of the smoothed time series with a temperature anomaly threshold ${\rm{\Delta }}{T}_{\mathrm{th}}$. As we are interested in the initial propagation front for each El Niño event, we additionally introduce an artificial refractory period of 11 months, which means threshold crossings occurring within 11 months after previous event are disregarded. Figure 16 depicts this procedure for a threshold value ${\rm{\Delta }}{T}_{\mathrm{th}}=2.5$ °C leading to the identification of three events in this example (vertical lines in figure 16). From this, we finally arrive at N = 140 discrete event series corresponding to the onsets of El Niño SST anomaly elevations at the different longitudinal locations along the dashed line in figure 15. The resulting spike trains are depicted in figure 17 on the left for four different threshold values ${\rm{\Delta }}{T}_{\mathrm{th}}=1.5$ °C–3.0 °C (step size of 0.5). The first thing to note is that this event detection procedure is able to correctly identify the El Niño occurrences as seen from the horizontal event structure in accordance with the El Niño years (labeled on the x-axis in the raster plots (a)–(d)) as well as the consistently high SPIKE-synchronization values C for all threshold values. With increasing temperature threshold, however, only the three strongest El Niño events are being captured (1983, 1998 and 2016, bold labels in figure 17) as seen from figures 17(c) and (d).

Zoom In Zoom Out Reset image size

To quantify the consistency of the propagation patterns, we perform a spike train order analysis on these event sequences. The results are also presented in figure 17. As a central observation, we found that for the three strong El Niño events (1983, 1998, 2016), there is a clear consistent propagation of the SST anomaly elevation at the equator from east to west as confirmed by a synfire indicator of ${F}_{{\rm{e}}-{\rm{w}}}=0.5$ for the east-to-west sorting (${\rm{\Delta }}{T}_{\mathrm{th}}=2.5$ °C, figure 17(c)). Furthermore, optimizing the synfire indicator by changing the sorting only leads to minor improvements (${F}_{{s}}=0.53$) and the resulting optimal order is very close to the original east-to-west sorting indicated in figure 17(g). Finally, this observation is also confirmed by the SPIKE-order based color coding in the raster plots (c) of the spike trains. For each of the three major El Niño events (1983, 1998, 2016) we see a clear tendency of leader (red) to follower (blue) going from east (top) to west (bottom) for all three threshold values. Further increase of the threshold temperature (figures 17(d) and (h)) leads to essentially the same observation, but parts of the longitudinal propagation might be missed as the SST anomaly maxima stay below threshold.

Moderate El Niño events (1992, 2003, 2010), on the other hand, seem to propagate in eastwards direction, as seen by opposite directed coloring for those years in figure 17(b) (${\rm{\Delta }}{T}_{\mathrm{th}}=2.0$). Consequently, the synfire indicator for the east-to-west sorting is very low, ${F}_{{\rm{e}}-{\rm{w}}}=0.12$, as some of the propagation events are moving in the opposite direction which reduces the synfire indicator. Therefore, it is also not possible to find a correct ordering, as indicated by a rather low optimized synfire indicator of ${F}_{{s}}=0.31$ in this case. This is even more pronounced for the smallest threshold temperature ${\rm{\Delta }}{T}_{\mathrm{th}}=1.5$, shown in figure 17(b). There, essentially all moderate El Niños are captured and most of those exhibit an eastwards propagation, opposed to the strong El Niños that propagate westwards. Therefore, the initial synfire indicator for the east-to-west sorting is negative, ${F}_{{\rm{e}}-{\rm{w}}}=-0.11$ and the optimal sorting shows the sign of a west-to-east propagation (figure 17(e)) with ${F}_{{s}}=0.38$.

To better understand the dependency of the spike train order analysis on the temperature threshold, we computed the synfire indicator for the initial east-to-west ordering as well as the value for optimal ordering ${F}_{{s}}$ for varying thresholds ${\rm{\Delta }}{T}_{\mathrm{th}}=1.5...3.5$ (figure 18). Additionally, figure 18 also shows the SPIKE-synchronization values for these parameters. One first observes that the SPIKE-synchronization remains rather constant over the whole range of threshold values, indicating a consistently good identification of the El Niño events. The synfire indicator for the original east-to-west sorting ${F}_{{\rm{e}}-{\rm{w}}}$ on the other hand shows a clear increase saturating in a plateau for a threshold value of around ${\rm{\Delta }}{T}_{\mathrm{th}}=2.5$ °C, where it is then also very close to the optimal value ${F}_{{s}}$. This is easily understood by the fact that at temperature thresholds above 2.5 °C, only the three strong El Niños are captured that show a consistent westwards propagation, while for values below also weaker El Niños that exhibit eastward propagation enter the analysis.

Zoom In Zoom Out Reset image size

In summary, we showed that the newly proposed spike train order method is also readily applicable to climate data and is able to identify propagation patterns in recurring events such as the El Niño. Our results indicate that close to the equator, strong El Niño events exhibit a clear westward propagation of an SST front, while for moderate El Niños we found indications for an eastward front propagation. Interestingly, previous results based on zonal current velocity within the El Niño 3.4 region found eastward propagation for the extreme El Niño events in 1983 and 1998 [24], while here we found west-ward propagation (figure 17). This discrepancy is most likely due to the different regions (El Niño 3.4 versus equatorial strip) and methodology (average warming/cooling rates versus front propagation).