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

Repetitive spatio-temporal propagation patterns are encountered in fields as wide-ranging as climatology, social communication and network science. In neuroscience, perfectly consistent repetitions of the same global propagation pattern are called a synfire pattern. For any recording of sequences of discrete events (in neuroscience terminology: sets of spike trains) the questions arise how closely it resembles such a synfire pattern and which are the spike trains that lead/follow. Here we address these questions and introduce an algorithm built on two new indicators, termed SPIKE-order and spike train order, that define the synfire indicator value, which allows to sort multiple spike trains from leader to follower and to quantify the consistency of the temporal leader-follower relationships for both the original and the optimized sorting. We demonstrate our new approach using artificially generated datasets before we apply it to analyze the consistency of propagation patterns in two real datasets from neuroscience (giant depolarized potentials in mice slices) and climatology (El Niño sea surface temperature recordings). The new algorithm is distinguished by conceptual and practical simplicity, low computational cost, as well as flexibility and universality.

In all of these fields recordings often exhibit well-defined patterns of spatio-temporal propagation where some prominent feature first appears at a specific location and then spreads to other areas until potentially becoming a global event.Such characteristic propagation patterns occur in phenomena such as avalanches (Lacroix et al. , 2012), tsunamis (Pelinovsky et al. , 2006), chemical waves and diffusion processes (Kuramoto , 2012), and epileptic Email address: thomas.kreuz@cnr.it(Thomas Kreuz).
In many cases spatio-temporal recordings can be represented as a two-dimensional plot where for each recording site the occurrence of certain discrete events (often obtained from threshold crossings in continuous data) are indicated by time markers.In neuroscience such a plot is known as a raster plot.A sequence of stereotypical neuronal action potentials (spikes, Rieke et al. (1996)) is a spike train and a set of spike trains exhibiting perfectly consistent repetitions of the same global propagation pattern is called a synfire pattern.In this paper we adapt this terminology and use all of these expressions not only in the context of neuronal spikes but also for any other kind of discrete events.However, note that our use of the term 'synfire pattern' differs slightly from the literature (see e.g.Abeles (2009)).Here we define a synfire pattern as a sequence of global events in which all neurons fire in consistent order and the interval between successive events is at least twice as large as the propagation time within an event.An example of a rasterplot with spike trains forming a perfect synfire pattern is shown in Fig. 1a.For any spike train set exhibiting propagation patterns the questions arise naturally whether these patterns show any consistency, i.e., to what extent do the spike trains resemble a synfire pattern, are there spike trains that consistently lead global events and are there other spike trains that invariably follow these leaders?Such questions about leader-follower dynamics have been specifically investigated not only in neuroscience (Pereda et al. , 2005), but also in fields as wide-ranging as, e.g., climatology (Boers et al. , 2014), social communication (Varni et al. , 2010), and human-robot interaction (Rahbar et al. , 2015).
In this study we introduce a framework consisting of two directional measures (SPIKE-Order and Spike Train Order ) that allows to define a value termed Synfire Indicator which quantifies the consistency of the leader-follower relationships in a rigorous and automated manner.This Synfire Indicator attains its maximal value of 1 for a perfect synfire pattern in which all neurons fire repeatedly in a consistent order from leader to follower (Fig. 1a).
The same framework also allows to sort multiple spike trains from leader to follower, as illustrated in Figs.1b and  1c.This is meant purely in the sense of temporal sequence.Whereas Fig. 1b shows an artificially created but rather realistic spike train set, in Fig. 1c the same spike trains have been sorted to become as close as possible to a synfire pattern.Now the spike trains that tend to fire first are on top whereas spike trains with predominantly trailing spikes are at the bottom.
We demonstrate the new approach using artificially generated datasets before we apply it to analyze the consistency of propagation patterns in two real datasets from very different fields of research, neuroscience and climatology.The neurophysiological dataset consists of neuronal activity recorded from mice brain slices.These recordings typi-cally exhibit a sequence of global events termed Giant Depolarized Potentials (GDPs) and one of the main questions we investigate is whether it is possible to identify neurons that consistently lead these events (potential hub neurons, see Bonifazi et al. (2009)).The climate dataset uses Optimum Interpolated Sea Surface Temperature (OISST) data provided by the National Oceanic and Atmospheric Administration (NOAA) to follow the El Niño phenomenon in the central pacific region Santoso et al. (2013) over the range of 35 years, from 1982 to 2016.We employ a threshold criterion to track the El Niño events and then quantify the consistency of the longitudinal movement of the propagation front.
The remainder of the article is organized as follows: In the Methods (Section 2) we first describe the coincidence detection (Section 2.1) and the symmetric measure SPIKE-Synchronization (Section 2.2).Subsequently, we introduce the new directionality approach consisting of the two measures SPIKE-Order and Spike Train Order (Section 2.3) as well as the Synfire Indicator (Section 2.4) before we discuss the use of SPIKE-Order surrogates to evaluate the statistical significance of the results in Section 2.5.The Results Section 3 consists of three Subsections detailing applications of the new approach to artificially generated datasets (Section 3.1), neurophysiological data (Section 3.2) and climate data from the El Niño phenomenon (Section 3.3).Conclusions are drawn in Section 4. Finally, both real datasets are described in the Appendix, the electrophysiological recordings in Appendix A, and the El Niño dataset) in Appendix B.

Measures
Analyzing leader-follower relationships in a spike train set requires a criterion that determines which spikes should be compared against each other.What is needed is a match maker, a method which pairs spikes in such a way that each spike is matched with at most one spike in each of the other spike trains.This match maker already exists.It is the adaptive coincidence detection first used as the fundamental ingredient for the bivariate measure event synchronization (Quian Quiroga et al. , 2002;Kreuz et al. , 2007a).
Event synchronization itself is symmetric and quantifies the overall level of synchrony from the number of quasi-simultaneous appearances of spikes.It was proposed along with an asymmetric measure termed delay asymmetry which evaluates the temporal order among matching spikes in the two spike trains.
However, unfortunately both event synchronization and delay asymmetry are defined for the bivariate case of two spike trains only, rely on sampled time profiles, and have a very non-intuitive normalization.For the symmetric variant we have already addressed these issues by proposing SPIKE-Synchronization (Kreuz et al. , 2015), a renormalized multivariate extension of event synchronization.
The two new measures SPIKE-Order and Spike Train Order proposed here improve and extend the asymmetric measure delay asymmetry in the same way.In particular, instead of just quantifying bivariate directionality they open up a completely new application, since they allow us to sort the spike trains according to the typical relative order of their spikes and to quantify the consistency of this order using the Synfire Indicator.
All four approaches (bivariate/multivariate, symmetric/asymmetric) are time-resolved as well as parameterand scale-free.Their calculation consists of two steps, adaptive coincidence detection followed by a combination of normalization and windowing.The first step, adaptive coincidence detection, is the same for all of these measures.

Adaptive Coincidence Detection
Most coincidence detectors rely on a coincidence window of fixed size τ (Kistler et al. , 1997;Naud et al. , 2011).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 Fig. 2a for an example), Quian Quiroga et al. (2002) proposed a more flexible coincidence detection.This coincidence detection is scale-and thus parameter-free since the minimum time lag τ (1,2) ij at which two spikes t (1) i and t (2) j of spike trains (1) and ( 2) are no longer considered to be synchronous is adapted to the local firing rates according to For some applications it might be appropriate to additionally introduce a maximum coincidence window τ 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.

SPIKE-Synchronization
In normalization and windowing SPIKE-Synchronization (Kreuz et al. , 2015) 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 (Quian Quiroga et al. , 2002), a more detailed comparison of the two measures can be found in Kreuz et al. (2015).
The main difference is that SPIKE-Synchronization (Kreuz et al. , 2015) results in a discrete, not a continuous, spike-timing based profile.The coincidence criterion is quantified by means of a coincidence indicator (2)

= min{t
(1) Thus any spike from the second spike train can at most be coincident with one spike from the first spike train.Small vertical lines mark the times right in the middle between two spikes, and a line is dashed when it does not mark the edge of a coincidence window.(c) In the same way a coincidence window τ from the second spike train.For two spikes to be coincident they both have to lie in each other's coincidence window which means that their absolute time difference has to be smaller than (2) j } (which is equivalent to the shorter definition found in Eq. 1).For the two spikes t (1) i and t (2) j on the left side this is the case, whereas the spikes on the right side are not coincident.
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 Fig. 2b, is the essential property which allows the adaptive coincidence detection to act as a match-maker for the subsequent application of SPIKE-Synchronization.Fig. 2c 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 Eq. 2 to all pairs of spike trains (n, m) with n, m = 1, ..., N and N denoting the number of spike trains: Here τ (n,m) ij is defined equivalent to Eq. 1. Subsequently, for each spike of every spike train a normalized coincidence counter 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: 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 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 Fig. 3a and 3b and one intermediate example (random distribution of spikes among spike trains) is shown in Fig. 3c.For a derivation of the expectation value for Poisson spike trains please refer to Mulansky et al. (2015).
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., 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 (n,m) i and the symmetric E (n,m) i , 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 The rasterplot shows 6 spike trains which emit spikes in nine reliable events.For the first two events spikes fire in order, for the next three events the order is random whereas for the last four events the order is inverted.In the last event there is one spike missing.Spike thickness decodes the SPIKE-Synchronization value C(t k ) (here almost constant), spike color the SPIKE-Order value D(t k ).(b,c) The SPIKE-Synchronization profile C(t k ) and its mirror profile (dashed black lines) act as envelope for both the SPIKE-Order profile D(t k ) (b, red) and the Spike Train Order profile E(t k ) (c, black).(b) The SPIKE-Order profile can not distinguish events with different firing order and by construction the average value is always D = 0. (c) On the other hand, in the Spike Train Order profile events with different firing order can clearly be distinguished.For the first two correctly ordered events the value 1 is obtained.The next three events exhibit random order and correspondingly the profile fluctuates rather wildly.Finally, the last four inversely ordered yield the value −1 except for the last event for which the absolute minimum value can not be obtained since one spike is missing.The average value, the Synfire Indicator F , is not 0 but negative which reflects the dominance of the inversely ordered events.

SPIKE-Order and Spike Train Order
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 symmet- where the index j is defined from the minimum in Eq. 2 as j = arg min j (t is obtained in an asymmetric manner as 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 (n,m) i = 0), but also in cases in which the times of the two coincident spikes are absolutely identical (t The multivariate profile D(t k ) obtained analogously to Eq. 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 Eqs.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 (Fig. 4a), it is not useful as an overall indicator of Spike Train Order (Fig. 4b).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).
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 Eqs.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: and 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 Eq. 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 Fig. 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 .

Synfire Indicator
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:  (n<m) , marked in black, is used to calculate the Synfire Indicator F , for both the unsorted spike trains (Fu, left) and the sorted spike trains (Fs, right).The thick black arrow in between the two matrices indicates the sorting process.
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.
It 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 Fig. 4c.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 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 (Fig. 5).
Hence if D (n,m) > 0 spike train n is leading m, while D (n,m) < 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) > 0 for n < m and D (n,m) < 0 for n > m.Therefore, we construct the overall SPIKE-Order as 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: The value is identical to the one of Eq. 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 ϕ(n).The overall Synfire Indicator for this permutation is then denoted as F ϕ .Accordingly, for the initial (unsorted) order of spike trains ϕ u the Synfire Indicator is denoted as The aim of the analysis is now to find the optimal (sorted) order ϕ s as the one resulting in the maximal overall Synfire Indicator F s = F ϕs : 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 ≈ 0.
The complexity of the problem to find the optimal Spike Train Order is similar to the well-known travelling salesman problem (Applegate et al. , 2011).For N spike trains there are N !permutations ϕ, so for large numbers of spike trains finding the optimal Spike Train Order ϕ s is a nontrivial 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 (Earl and Deem , 2005) or simulated annealing (Dowsland and Thompson , 2012) 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 ∆F = −2D (k,k+1) .All moves with positive ∆F are accepted while the likelihood of accepting moves with negative ∆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 Fig. 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 (Eq.14, Fig. 6c) and the average value of the Spike Train Order profile E(t k ) (Eq.11, Fig. 6d).Finally, the sorted spike trains in Fig. 6e 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 ).
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, For the first two events the synfire-order of the original is destroyed in the surrogates whereas for the next two events both sequences are equally unordered.(b) Histogram for 19 surrogates.Thick lines denote mean and standard deviation.Since the value for the original dataset (black) is not maximum, the optimally sorted spike trains do not exhibit a statistically significant 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.

Statistical significance
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 (Fig. 7a) 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 Fig. 7a 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 Eq. 1 and change the value of SPIKE-Synchronization.This in turn would make the desired evaluation of pure spike order effects difficult.
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 σ: Results of the significance analysis for our standard example are shown in the histogram in Fig. 7b.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.

Results
In the following we apply our new algorithm to artificially generated datasets (Section 3.1), neurophysiological data (Section 3.2) and, finally, climate data from the El Niño phenomenon (Section 3.3).For the Figures we use the same full layout introduced in Fig. 6 to which we add the significance analysis of Fig. 7b.

Application to artificially generated data
We start with examples covering the two extreme cases of a perfect synfire pattern and a completely random spike train set.First, in Fig. 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 Fig. 8. SPIKE-and Spike Train Order analysis for a perfect inverse synfire pattern.The plot follows the layout of Fig. 6 with the histogram of the surrogate test (see Fig. 7b) for statistical significance added as subplot f.For the unsorted spike trains a minimal Synfire Indicator of Fu = −1 is obtained, while sorting results in the maximum value of Fs = 1.According to the surrogate test the statistical significance of the result is very high.
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 (Fig. 8f) shows that the statistical significance of the original Synfire Indicator is very high.The other extreme case is Poisson spike trains (Fig. 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.
The third example in Fig. 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 Fig. 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 Fig. 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 Fig. 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 preserv- Fig. 10.SPIKE-and Spike Train Order analysis for Poisson spike train interspersed with spike trains that contain random spikes but also a perfect inverse synfire pattern.The order contained within the synfire pattern spike train is distinct enough to make the Synfire Indicator for the sorted spike trains statistically significant.
ing the spike numbers in the surrogates it explicitly takes the statistics of each dataset into account.

Application to neuroscience
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 Ben-Ari et al. (1989)).In previous studies, GDPs have been used to investigate the topology of networks and the role of hub cells (Bonifazi et al. , 2009) as well as to reveal the deterministic and stochastic processes underlying spontaneous, synchronous network bursts (Takano et al. , 2012).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 Fig. 11 includes 13 GDPs to their firing rate starting with the most sparse spike trains.The messy color-patterns reveal that this is completely uncorrelated to the spike order within the GDPs.(f) After sorting, there is a fairly consistent transition from spike trains with predominantly leading spikes (red) in the GDPs to spike trains with predominantly following spikes (blue).
over a bit more than 6 minutes.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 (Fig. 11a).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 Fig. 11b 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 Fig.12. SPIKE-Order for the real data already analyzed in Fig. 11 but this time the analysis of SPIKE-Order was restricted to spikes with a SPIKE-Synchronization value of at least 0.7.This simple thresholding allows to focus the analysis on the reliable events and to disregard all spikes between the events (these are not colored and thus remain black).This results in an increase of the overall value of SPIKE-Order from 0.284 to 0.438.
be seen by comparing the color bars on the left of Fig. 11a and Fig. 11e.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 (Fig. 11d).Finally, the surrogate analysis (Fig. 11f) shows this result to be highly significant.However, the spiking in Fig. 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 SPIKEorder 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 Fig. 12 for the same dataset already used in Fig. 11.Focusing the analysis on the reliable GDPs leads to an increase of the Synfire Indicator from 0.284 to 0.438.
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 Fig. 13 where we have color-coded the optimized Spike Train Order obtained in Fig. 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.
In Fig. 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.
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 Fig. 11) and 0.92 after (as in Fig. 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, Fig. 14.SPIKE-Order for a second dataset for which again the analysis of SPIKE-Order was restricted to spikes with a SPIKE-Synchronization value of at least 0.7.In this case the focus on the GDPs let the Synfire Indicator increase from 0.296 (result not shown) to 0.389.In addition, here we also calculated one average Spike Train Order value per GDP (green points) which illustrates once more the time-resolved nature of the method.
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 Takano et al. (2012).

Application to climate data
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.Malik et al. (2012); Boers et al. (2014).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 (McPhaden et al. , 2006) and therefore this phenomenon has been studied extensively in terms of both intensity (e.g.Rasmusson et al. (1982); Yeh et al. (2009)) and frequency (e.g.Qian et al. (2011); An and Wang (2000)).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 (Santoso et al. , 2013;Antico and Barros , 2016).The region used here to analyze the SST is depicted in Fig. 15 (dashed line), and an examplary smoothed time profile for the center of the analyzed region (215 • E) is shown in Fig. 16.Further details on the region, the dataset and data preprocessing (smoothing) is outlined in Appendix B.
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 ∆T 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.Fig. 16 depicts this procedure for a threshold value ∆T th = 2.5 • C leading to the identification of three events in this example (vertical lines in Fig. 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 Fig. 15.The resulting spike trains are depicted in Fig. 17 on the left for four different threshold values ∆T 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 ad) 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 Fig. 17) as seen from Fig. 17c and Fig. 17d.
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 Fig. 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 e-w = 0.5 for the east-to-west sorting (∆T th = 2.5 • C, Fig. Fig. 17c).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 Fig. 17g.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 (Fig. 17d 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 Fig. 17b (∆T th = 2.0).Consequently, the Synfire Indicator for the east-to-west sorting is very low, F e-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 ∆T th = 1.5, shown in Fig. 17b.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 e-w = −0.11and the optimal sorting shows the sign of a west-to-east propagation (Fig. 17e) with F s = 0.38.
To better understand the dependency of the Spike Train Order analysis on the temperature threshold, we computed  -d) show the event raster plots with with SPIKE-Order value in color coding using an eastto-west ordering.The years shown on the x-axis (representing January 1st of the indicated year) denote the moderate and strong (bold face) El Niño events.On the right (e-h), the optimal leader-follower order compared to an initial east-to-west sorting is shown (circles: spike trains with at least 2 spikes, triangles: spike trains with less than two spikes).Additionally, the SPIKE-Synchronization C as well as the Synfire Indicator for the original (east-to-west) Fe-w, and the optimal sorting Fs are reported.the Synfire Indicator for the initial east-to-west ordering as well as the value for optimal ordering F s for varying thresholds ∆T th = 1.5 . . .3.5 (Fig. 18).Additionally, Fig. 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 e-w on the other hand shows a clear increase saturating in a plateau for a threshold value of around ∆T 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.
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 (Santoso et al. , 2013), while here we found west-ward propagation (Fig. 17).This discrepancy is most likely due to the different regions (El Niño 3.4 vs equatorial strip) and methodology (average warming/cooling rates vs front propagation).

Discussion
Over the last years a wide variety of measures to quantify the synchrony between spike trains have been introduced.Three recent proposals, ISI-distance (Kreuz et al. , 2007b(Kreuz et al. , , 2009)), SPIKE-distance (Kreuz et al. , 2011(Kreuz et al. , , 2013)), and SPIKE-Synchronization (Kreuz et al. , 2015;Mulansky et al. , 2015), share the desirable property of being time-resolved and parameter-free (time-scale independent).However, their bivariate versions are symmetric and in consequence their multivariate versions are invariant to changes in the order of spike trains.None of these measures is designed to provide information about the directionality of the propagation patterns.
In the present study we address this issue.First we use an adaptive coincidence detection as match maker in order to identify pairs of coincident spikes.Then we define two measures, the asymmetric SPIKE-Order D and the symmetric Spike Train Order E, which are particularly useful in a bivariate representation (pairwise matrix) and as a time-resolved multivariate profile, respectively.From these two measures we can derive the Synfire Indicator F , a condensed scalar value that quantifies the overall consistency of the spatio-temporal propagation patterns in a rigorous and automated way.Its maximization allows to sort multiple spike trains from leader to follower.This is meant purely in the sense of temporal sequence of the spikes.The question asked is: For which spike trains do spikes tend to occur first and for which do they tend to occur last?We use simulated annealing to search among all permutations of spike trains for the sorting that resembles as closely as possible a synfire pattern, a perfectly consistent repetition of the same global propagation pattern.In a final step we evaluate the statistical significance of the optimized permutation using a set of carefully constructed spike train surrogates.
We first illustrate the merits of our new approach using artificially generated datasets and then apply it to real datasets from two very different fields of research, neuroscience and climatology.In the neurophysiological dataset we analyze Giant Depolarized Potentials (GDPs) recorded from mice brain slices in order to search for potential hub neurons (Bonifazi et al. , 2009), whereas in the climate application we quantify the consistency of the longitudinal movement of the propagation front of El Niño events (Santoso et al. , 2013).
The new algorithm is conceptually simple, of low computational cost and comes with an intuitive and straightforward visualization, including a color-coded rasterplot.It substantially improves on all the bivariate functionalities of its predecessor directional measure delay asymmetry (no need for sampled profiles, more intuitive normalization etc.) and could thus also be used in the context of the pairwise matrices, both normalized or cumulative, used in complex network theory (Malik et al. , 2010;Boers et al. , 2014).However, one of the main advantages of the new algorithm is its multivariate nature which opens up completely new kinds of application such as spike train sorting.
One important advantage that the method shares with other techniques of spike train analysis is the high flexibility in the definition of events.For example, when looking at the synchronization of neuronal bursts instead of individual spikes one can define the events as the onset, the center or the offset of the activity (e.g., the first, the middle or the last spike of each burst).In cases in which a burst of spikes is considered to be equivalent to a single spike one could introduce some kind of meta-events and then look at coincidences between these meta-events.
The application of our measures is also not restricted to truly discrete data.Continuously sampled data can be re-duced to a spike train where the only information maintained is the timing of the individual events.Often these event times are obtained in a manner similar to how the neuronal spike times are extracted from recordings of neuronal membrane potentials (usually done via some kind of thresholding).Examples of sampled data to which measures of spike train synchrony have been applied include EEG data (Quian Quiroga et al. , 2002;Kreuz et al. , 2013;Rosário et al. , 2015) and, outside of neuroscience, stock market velocity (Zhou and Sornette , 2003) -and rainfall events (Malik et al. , 2010).
The algorithm is particularly suited for datasets with a high value of SPIKE-Synchronization.According to the coincidence criterion (Eq. 1) these are spike trains that include sequences of global events for which the interval between successive events is at least twice as large as the propagation time within an event.For these datasets the Synfire Indicator evaluates the consistency of the order within these well separated global events.The universality of the phenomenon, repetitive propagation patterns, makes our new algorithm applicable in a wide array of fields such as medical sciences, seismology, oceanography, meteorology or climatology.For example, the duration of an epileptic seizure is typically much shorter than the interval between two successive seizures.Also the time it takes a storm front to cross a specific region is typically much smaller than the time to the next storm.Many other repetitive propagation phenomena exhibit similar ratios of characteristic time scales.
In order to understand the scope of our proposed algorithm it is important to understand what it is not designed to achieve.The method deals purely with relative order, it does not consider the length of absolute delays.Moreover, while the instantaneous coincidence criterion makes the method time-scale independent, parameter-free and easy to use, it also renders it insensitive to patterns involving spikes that are not immediately adjacent.Many other, typically more complicated, methods have been designed to characterize the detailed spatio-temporal structure in large neuronal networks (Schneider et al. , 2006;Gansel and Singer, 2012) or to detect hierarchically structured spike-train communities (Billeh et al. , 2014;Russo and Durstewitz, 2017).The method is also not designed to detect neuronal synfire chains (in the strict sense of the word) in massively parallel data.For this task other statistical methods based on some forms of pattern detection have been developed (Schrader et al. , 2008;Gerstein et al. , 2012).
Another caveat concerns causality.While a significant value of the Synfire Indicator F s in our algorithm clearly shows the presence of a preferred temporal order of some signals with respect to others, it does not necessarily prove a driver-responder relationship.There are other methods that have been developed for this kind of system dynamics analysis (e.g., Andrzejak and Kreuz (2011)).But even for such methods causality is always a strong claim.In fact, the two signals might be driven by a common hidden source and a consistent leader (follower) could just indicate a drive with a smaller (larger) delay.Similarly, internal delay loops in one of the two systems can also fool the interpretation.
There are a number of possible directions for future research, both from a methods and from a data point of view.Regarding the algorithm, for the coincidence detection it would be straightforward to limit the range of allowed time lags by incorporating information about the expected speed of propagation (Malik et al. , 2010).One could introduce a minimum time lag in order to ensure causality and/or limit the maximally allowed time lag in order to focus on meaningful propagation of activity.In principle the range of allowed time lags could even be selected individually for each pair of spike trains depending on the known properties of the connectivity between the respective two neurons.Importantly, even with such type of time lag restrictions in place, it has still to be guaranteed that each spike can be part of at most one coincidence.
A follow-up task for our neurophysiological data would be to investigate to what extent the neurons that are identified as leading by our analysis are identical to the so-called hub neurons (Bonifazi et al. , 2009), i.e. neurons with a much higher than average degree of connectivity within the network.Regarding the El Niño analysis, the difference in propagation directions observed for the strong and weak El Niño events remains an open question.Furthermore, expanding the analysis to wider regions, i.e. 5 • S . . . 5 • N and verifying the consistency of the observed propagation patterns is an interesting topic for future research.Finally, the relation to results based on average warming/cooling rates (Santoso et al. , 2013) requires further investigations.
The algorithm will be readily applicable for everyone since, together with the existing symmetric measures ISIdistance, SPIKE-distance, and SPIKE-Synchronization, SPIKE-Order is implemented in in three publicly available software packages, the Matlab-based graphical user interface SPIKY1 (Kreuz et al. , 2015), cSPIKE2 (Matlab command line with MEX-files), and the Python library PySpike3 (Mulansky et al. , 2016).

Appendix A. Neurophysiological dataset
The neurophysiological data analyzed in Section 3.2 were recorded via fast multicellular calcium imaging in acute CA3 hippocampal slices from juvenile mice.The CA3 region has a strong recurrent excitatory connectivity (Amaral and Witter, 1995).This distinct feature is suggested to be crucial for memory encoding and pattern completion and thus memory retrieval (Leutgeb and Leutgeb, 2007).During memory retrieval in rodents, population bursts of the CA3 lead to high frequency stimulation of the efferent regions, so called sharp wave ripples (Chrobak and Buzsáki, 1996).In the juvenile hippocampus, due to a higher chloride reversal potential in the CA3 pyramidal cells, the GABAergic system is excitatory (Ben-Ari et al. , 1989).GABAergic interneurons have been shown to serve as so called hub neurons that trigger the GDPs (Bonifazi et al. , 2009).
The recordings were performed by the group of Heinz Beck at the Department of Epileptology, University of Bonn, Germany, prior to and independently from the design of this study.Transversal acute brain slices (300µm thick) were prepared from 5 to 10-day-old (P5-P10) C57BL/6 mice (Charles River, n = 19 slices).Slice preparation, calcium imaging and data analysis were performed as previously described in (Hedrich et al. , 2014).For AMloading of brain slices with OGB1-AM we used a protocol modified from (Crépel et al. , 2007).Multicellular calcium imaging was done using a homemade single planar illumination microscope (SPIM) modified from (Holekamp et al. , 2008).Movies were recorded at a frame rate of 200 Hz over a minimal length of 5 min up to 30 min to record a sufficient amount of spontaneous activity.Time points of cell activity from the imaging data were defined as the onsets of Ca 2 events in fluorescence traces of all individual cells using the maximum of the second derivative of each event (Henze et al. , 2000).
In order to test the Spike Train Order algorithm, datasets were chosen that exhibited at least three global GDPs during the recording (n = 5).For one dataset the surrogate analysis described in Section 2.5 proved to be unfeasible due to its excessive density of spontaneous activity.Therefore this dataset was discarded from further analysis, so the final number of datasets analyzed was n = 4.

Appendix B. El Niño dataset
The data analyzed in Section 3.3 describe one of the most well-known global climate effects, the so called El Niño, the warm phase of the El Niño Southern Oscillation (ENSO).It is usually characterized by an increased sea surface temperature (SST) in the central and eastern tropical Pacific Ocean typically during the months September -February (Trenberth , 1997).According to the U.S. National Oceanic and Atmospheric Administration (NOAA), an El Niño event is identified by an increased three-month moving average SST by 0.5 • C over at least six months.El Niños can last from nine months up to two years and typically occur in irregular intervals of two to seven years.
The analysis presented in Section 3.3 uses the Optimum Interpolated Sea Surface Temperature (OISST) data provided by NOAA.These data result from a high-resolution blended analysis (spatial resolution of 0.25 • ) of daily SST and ice constructed by combining observations from different platforms (satellites, ships, buoys) on a regular global grid with a time range from 1981 to 2016 (Reynolds et al. , 2007).The daily SST data is centralized by the long-term daily mean resulting in daily SST anomaly data (deviations from long-term mean) that form the basis of this analysis.
The area used to define El Niño events, the El Niño 3.4 region, stretches from 5 • S . . . 5 • N in latitude and 190 • E . . .240 • E in longitude, as shown in Fig. 15.However, for the latitudinal direction we here focused only on a small central strip around the equator, i.e. from 0.5 • S . . .0.5 • N over a longitudinal stretch consistent with the El Niño region, i.e. from 180 • E . . .250 • E (dashed line in Fig. 15).Note, that this focus on the small region around the equator is the main difference to most previous works that studied the propagation of SST anomalies averaged over a much larger latitudinal extent, e.g. 5 • S . . . 5 • N in Santoso et al. (2013); Antico and Barros (2016).In latitudinal direction, we average over the whole strip (four grid points), while in longitudinal direction an averaging over two grid points is performed resulting in 0.5 • of spatial resolution and a total of N = 140 time series of daily SST anomalies.We disregard short-term fluctuation by applying a Gaussian smoothing with width σ = 14 days.In Fig. 16 we show the time series resulting from this procedure exemplarily for the center of the observed region, 215 • E.
For the Sea Surface Temperature data we use a threshold criterion to identify the moving high temperature fronts of the El Niño events.The threshold determines the signal-tonoise ratio of the data: for small values even weaker events have an effect on the result, whereas for higher values the analysis is focused on the strongest El Niño events only.Due to the variability of the propagation patterns it might happen that in successive years the threshold is surpassed in different regions.However, since the aim of the analysis is to look at the propagation of individual (seasonal) fronts we suppress coincidences between threshold crossings from different years.We still use the adaptive coincidence detection from Eq. 1 but define a maximum coincidence window τ max which in this case is set to 9 months.

Fig. 1 .
Fig. 1.Motivation for SPIKE-order and Spike Train Order.(a) Perfect synfire pattern.(b) Unsorted set of spike trains.(c) The same spike trains as in B but now sorted from leader to follower.

Fig. 2 .
Fig. 2. (a) This example demonstrates the usefulness of an adaptive coincidence detection.Depending on context the same two spikes (left) can appear as coincident (right, top) or as non-coincident (right, bottom).(b) Illustration of the adaptive coincidence detection.For clarity spikes and their coincidence windows are shown alternatively in bright and dark color.The first step assigns to each spike t (1) i of the first spike train a potential coincidence window which does not overlap with any other coincidence window: τ (1) i

Fig. 3 .
Fig. 3. SPIKE-Synchronization.Note that the profile C(t k ) is defined only at the times of the spikes but a better visualization is achieved by connecting the individual dots.By construction the pooled spike train of these examples is identical consisting of 10 evenly spaced bursts.The only difference is the distribution of the spikes among the individual spike trains which varies from maximum to minimum via intermediate synchrony.SPIKE-Synchronization correctly indicates these changes.(a) Maximum reliability results in the value one over the whole time interval.Each spike train contains one spike per firing event.(b) Synfire pattern of bursts resulting in minimum reliability corresponding to the value zero for the whole time interval.(c) A random distribution of spikes among spike trains yields intermediate values.
SPIKE-Synchronization assigns to each spike of a given spike train pair a bivariate coincidence indicator.These coincidence indicators C (n,m) i

Fig. 4 .
Fig.4.SPIKE-Order profile D(t k ) and Spike Train Order profile E(t k ) for an artificially created example dataset.(a) The rasterplot shows 6 spike trains which emit spikes in nine reliable events.For the first two events spikes fire in order, for the next three events the order is random whereas for the last four events the order is inverted.In the last event there is one spike missing.Spike thickness decodes the SPIKE-Synchronization value C(t k ) (here almost constant), spike color the SPIKE-Order value D(t k ).(b,c) The SPIKE-Synchronization profile C(t k ) and its mirror profile (dashed black lines) act as envelope for both the SPIKE-Order profile D(t k ) (b, red) and the Spike Train Order profile E(t k ) (c, black).(b) The SPIKE-Order profile can not distinguish events with different firing order and by construction the average value is always D = 0. (c) On the other hand, in the Spike Train Order profile events with different firing order can clearly be distinguished.For the first two correctly ordered events the value 1 is obtained.The next three events exhibit random order and correspondingly the profile fluctuates rather wildly.Finally, the last four inversely ordered yield the value −1 except for the last event for which the absolute minimum value can not be obtained since one spike is missing.The average value, the Synfire Indicator F , is not 0 but negative which reflects the dominance of the inversely ordered events.

Fig. 5 .
Fig. 5. Pairwise cumulative SPIKE-Order matrix D before (left) and after (right) sorting for the example dataset from Fig. 4. The upper triangular matrix D(n<m) , marked in black, is used to calculate the Synfire Indicator F , for both the unsorted spike trains (Fu, left) and the sorted spike trains (Fs, right).The thick black arrow in between the two matrices indicates the sorting process.

Fig. 6 .
Fig. 6.Complete illustration of SPIKE-Order using our example dataset from Fig. 4. (a) Unsorted spike trains with the spikes color-coded according to the value of the SPIKE-Order D(t k ).(b) Spike Train Order profile E(t k ).The Synfire Indicator Fu for the unsorted spike trains is slightly negative.(c) Pairwise SPIKE-Order matrix D before and after sorting.The optimal order maximizes the upper triangular matrix.(d) Spike Train Order profile E(t k ) and its average values, the Synfire Indicator Fs for the sorted spike trains.(e) Sorted spike trains.

Fig. 7 .
Fig. 7. Statistical significance: Surrogate analysis for the example dataset from Fig. 4. (a) Spike order patterns for original (black) and one randomized surrogate (red).For clarity only the first four events are shown.For the first two events the synfire-order of the original is destroyed in the surrogates whereas for the next two events both sequences are equally unordered.(b) Histogram for 19 surrogates.Thick lines denote mean and standard deviation.Since the value for the original dataset (black) is not maximum, the optimally sorted spike trains do not exhibit a statistically significant synfire pattern.

Fig. 9 .
Fig.9.SPIKE-and Spike Train Order analysis for 20 Poisson spike trains.Since the number of spike trains is too large to label the spike trains in the top and in the bottom subplot with numbers we use color coding at the left side to label them.Both before and after sorting the Synfire Indicator is very close to zero.The surrogate analysis reveals the result to be non-significant.

Fig. 11 .
Fig. 11.SPIKE-Order for real data recorded in an acute hippocampal slice from a juvenile mouse.Note how the color-coding of the spikes according to their SPIKE-order D helps to overcome the low temporal resolution of the Figure and to resolve the spike order within the GDPs.(a) Initially the spike trains are sorted according to their firing rate starting with the most sparse spike trains.The messy color-patterns reveal that this is completely uncorrelated to the spike order within the GDPs.(f) After sorting, there is a fairly consistent transition from spike trains with predominantly leading spikes (red) in the GDPs to spike trains with predominantly following spikes (blue).

Fig. 13 .
Fig. 13.Projection of the optimized Spike Train Order on the 2D-photo of the hippocampal slice.The Regions of Interest (ROIs) which denote filled and identified cells in the CA3 region are colorcoded from leader (index 1, red) to follower (index 163, blue) using the optimized Spike Train Order of Fig. 12.The very first leader (lower right) and the very last follower (upper left) are marked by filled contours.
Fig. 15.El Niño 3.4 region (rectangle) and the small strip around the equator used for the analysis in this work (dashed line).Color coding represents daily SST anomaly (in • C) on November 11, 2015.
Fig.18.SPIKE-Synchronization C and Synfire Indicator for the original sorting (east-to-west) Fe-w and the optimal sorting Fs in dependence of the temperature threshold ∆T th .
• S . . .0.5 • N, 170 • E . . .250 • E) smoothed in time using a Gaussian kernel with width of 2 weeks.Events are defined as initial threshold crossings of the smoothed temperature time series (see text).From top to bottom four different threshold values have been used: ∆T th = 1.5 • C . . .3.0 • C. The left graphs (a