Measures of spike train synchrony and directionality
aa r X i v : . [ phy s i c s . d a t a - a n ] J a n Measures of spike train synchrony anddirectionality
Eero Satuvuori and Irene Malvestio and Thomas Kreuz
Abstract
Measures of spike train synchrony have become important tools in bothexperimental and theoretical neuroscience. Three time-resolved measures called theISI-distance, the SPIKE-distance, and SPIKE-synchronization have already beensuccessfully applied in many different contexts. These measures are time scale in-dependent, since they consider all time scales as equally important. However, in realdata one is typically less interested in the smallest time scales and a more adaptiveapproach is needed. Therefore, in the first part of this Chapter we describe recentlyintroduced generalizations of the three measures, that gradually disregard differ-ences in smaller time-scales. Besides similarity, another very relevant property ofspike trains is the temporal order of spikes. In the second part of this chapter weaddress this property and describe a very recently proposed algorithm, which quan-tifies the directionality within a set of spike train. This multivariate approach sortsmultiple spike trains from leader to follower and quantifies the consistency of thepropagation patterns. Finally, all measures described in this chapter are freely avail-able for download.
Eero SatuvuoriDepartment of Physics and Astronomy, University of Florence, Sesto Fiorentino, ItalyFaculty of Behavioural and Movement Sciences ,Vrije Universiteit Amsterdam, Netherlandse-mail: [email protected]
Irene MalvestioDepartment of Information and Communication Technologies, Universitat Pompeu Fabra,Barcelona, SpainDepartment of Physics and Astronomy, University of Florence, Sesto Fiorentino, Italye-mail: [email protected]
Thomas KreuzInstitute for Complex Systems, CNR, Sesto Fiorentino, Italye-mail: [email protected]
The brain can be considered as a huge network of spiking neurons. It is typicallyassumed that only the spikes, and not the shape of the action potential nor thebackground activity, convey the information processed within this network [11].Sequences of consecutive spikes are called spike trains. Measures of spike trainsynchrony are estimators of the similarity between two or more spike trains, whichare important tools for many applications in neuroscience. Among others, they al-low to test the performance of neuronal models [10], they can be used to quantifythe reliability of neuronal responses upon repeated presentations of a stimulus [18],and they help in the understanding of neural networks and neural coding [27].Over the years many different methods have been developed in order to quantifyspike train synchrony. They can be divided in two classes: time-scale dependent andtime-scale independent methods. The two most known time-scale dependent meth-ods are the Victor-Purpura distance [28] and the van Rossum distance [25]. Theydescribe spike train (dis)similarity based on a user-given time-scale to which themeasures are mainly sensitive to. Time scale independent methods have been devel-oped more recently. In particular, the ISI-distance [15], the SPIKE-distance [12, 13]and SPIKE-synchronization [16] are parameter-free distances, with the capability ofdiscerning similarity across different spatial scales. All of these measures are time-resolved, so they are able to analyze the time dependence of spike train similarity.One problematic aspect of time-scale independent methods is that they con-sider all time-scales as equally important. However, in real data one typically isnot interested in the very small time scales. Especially in the presence of bursts(multiple spikes emitted in rapid succession), a more adaptive approach that gradu-ally disregards differences in smaller time-scales is needed. Thus, in the first partof this chapter we describe the recently developed adaptive extensions of thesethree parameter-free distances: A-ISI-distance, A-SPIKE-distance and A-SPIKE-synchronization [23].All of these similarity measures are symmetric and in consequence invariant tochanges in the order of spike trains. However, often information about directionalityis needed, in particular in the study of propagation phenomena. For example, inepilepsy studies, the analysis of the varying similarity patterns of simultaneouslyrecorded ensembles of neurons can lead to a better understanding of the mechanismsof seizure generation, propagation, and termination [24, 4].In the second part of this chapter we address the question: Which are the neuronsthat tend to fire first, and which are the ones that tend to fire last? We present SPIKE-Order [17], a recently developed algorithm which is able to discern propagationpattern in neuronal data. It is a multivariate approach which allows to sort multiplespike trains from leader to follower and to quantify the consistency of the temporalleader-follower relationships. We close this chapter by describing some applicationsof the methods presented. easures of spike train synchrony and directionality 3
Two of the most well known spike train distances, the Victor-Purpura [28] and thevan Rossum distance [25], are time-scale dependent. One drawback of these meth-ods is the fixed time-scale, since it sets a boundary between rate and time coding forthe whole recording. In the presence of bursts, where multiple spikes are emitted inrapid succession, there are usually many time-scales in the data and this is difficultto detect when using a measure that is sensitive to only one time-scale at a time [6].The problem of having to choose one time-scale has been eliminated in thetime-scale independent ISI-distance [15], SPIKE-distance [12, 13] and SPIKE-synchronization [16], since these methods always adapt to the local firing rate. TheISI-distance and the SPIKE-distance are time resolved, time-scale free measures ofdissimilarity between two or more spike trains. The ISI-distance is a measure ofrate dissimilarity. It uses the interspike intervals (ISIs) to estimate local firing rateof spike trains and measures time-resolved differences between them. The SPIKE-distance, on the other hand, compares spike time accuracy between the spike trainsand uses local firing rates to adapt to the time-scale. SPIKE-synchronization is alsotime-scale free and is a discrete time resolved measure of similarity based on ISIderived coincidence windows that determine if two spikes in a spike train set arecoincident or not.The ISI-distance, SPIKE-distance, and SPIKE-synchronization are looking at alltime-scales at the same time. However, in real data not all time-scales are equallyimportant, and this can lead to spuriously high values of dissimilarity when lookingonly at the local information. Many sequences of discrete events contain differ-ent time-scales. For example, in neuronal recordings besides regular spiking oneoften finds bursts, i.e., rapid successions of many spikes. The A-ISI-distance, A-SPIKE-distance and A-SPIKE-synchronization [23] are generalized versions of pre-viously published methods the ISI-distance [14], SPIKE-distance [12] and SPIKE-synchronization [16]. The generalized measures also contain a notion of global con-text that discriminates between relative importance of differences in the global scale.This is done by means of a normalization based on a minimum relevant time-scale(MRTS). They start to gradually ignore differences between spike trains for inter-spike intervals (ISIs) that are smaller than the MRTS. The generalization providedby the MRTS is implemented with the threshold parameter thr , which is then ap-plied in a different way to each of the measures. The threshold is used to determineif a difference between the spike trains should be assessed in a local or in a globalcontext. This threshold is used for all three measures, but the way it is applied varies.The extended methods fall back to the original definitions when thr = N , indices ofspike trains by n and m , spike indexes by i and j and the number of spikes in spiketrain n by M n . The spike times of spike train n are denoted by { t ( n ) i } with i = . . . M n . Eero Satuvuori and Irene Malvestio and Thomas Kreuz mn (a) S p i k e t r a i n s t t (n)i t (n)i+1 t (m)j t (m)j+1 x ISI(n) (t)x
ISI(m) (t) mn (b) S p i k e t r a i n s t x (n)P (t) x (n)F (t)x (m)P (t) x (m)F (t) ∆ t i(n) ∆ t i+1(n) ∆ t j(m) ∆ t j+1(m) mn (c) t i(n) τ i(n) t j(m) τ j(m) t i(n) τ i(n) t j(m) τ j(m) Time S p i k e t r a i n s Fig. 1
Schematic drawing for all three measures. (a) Illustration of the variables used to define the
A-ISI-distance . All measures use the instantaneous interspike interval x ( n ) ISI ( t ) to adapt to the localfiring rate. (b) Additional variables used for the A-SPIKE-distance . (c) Coincidence criterion forthe
A-SPIKE-synchronization . The coincidence window of each spike is derived from the twosurrounding interspike intervals. For simplicity the thr = The A-ISI-distance [23] measures the instantaneous rate difference between spiketrains (see Fig. 1a). It relies on a time-resolved profile, meaning that in a first stepa dissimilarity value is assigned to each time instant. To obtain this profile, we firstassign to each time instant t the times of the previous spike and the following spike easures of spike train synchrony and directionality 5 S p i k e t r a i n s (a) I A ( t ) D IA =0.307 (b) S A ( t ) D SA =0.193 (c) Time [a.u.] C A ( t k ) S CA =0.480 (d) Fig. 2
Profiles of
A-ISI-distance (a),
A-SPIKE-distance (b) and
A-SPIKE-synchronization (c)for an artificial example dataset of 50 spike trains with population events with different jitters anddecreasing noise over time. t ( n ) P ( t ) = max { t ( n ) i | t ( n ) i ≤ t } for t ( n ) t t ( n ) M n (1) t ( n ) F ( t ) = min { t ( n ) i | t ( n ) i > t } for t ( n ) t t ( n ) M n . (2)From this for each spike train n an instantaneous ISI can be calculated as x ( n ) ISI ( t ) = t ( n ) F ( t ) − t ( n ) P ( t ) . (3)The A-ISI-profile is defined as a normalized instantaneous ratio in ISIs: I An , m ( t ) = | x ( n ) ISI ( t ) − x ( m ) ISI ( t ) | max { x ( n ) ISI ( t ) , x ( m ) ISI ( t ) , thr } . (4)For the A-ISI-distance the MRTS is defined so that when the ISI of both spiketrains are smaller than a threshold value thr , the threshold value is used instead. Themultivariate A-ISI-profile is obtained by averaging over all pairwise A-ISI-profiles: I A ( t ) = N ( N − ) / N − ∑ n = N ∑ m = n + I An , m ( t ) . (5)This is a non-continuous piecewise constant profile and integrating over time givesthe A-ISI-distance: Eero Satuvuori and Irene Malvestio and Thomas Kreuz D AI = t e − t s Z t e t s I A ( t ) dt . (6)Where t s and t e are the start and end times of the recording respectively. If thr is setto zero, the method falls back to the ISI-distance [14].Fig. 2a shows an artificial spike train dataset together with the correspondingA-ISI-profile in Fig. 2b. The A-ISI-profile for the example dataset shows high dis-similarity for the left side of the raster plot, where noise is high. When the noise isdecreased and rates become more similar in the right side, the dissimilarity profilegoes down. The overall ISI-distance is the mean value of the profile. The A-SPIKE-distance [23] measures the accuracy of spike times between spiketrains relative to local firing rates (see Fig. 1b). In order to assess the accuracy ofspike events, each spike is assigned a distance to its nearest neighbour in the otherspike train: ∆ t ( n ) i = min j ( | t ( n ) i − t ( m ) j | ) . (7)The distances are interpolated between spikes using for all times t the time differ-ences to the previous and to the following spikes x ( n ) P ( t ) and x ( n ) F ( t ) : x ( n ) P ( t ) = t − t ( n ) i for t ( n ) i t t ( n ) i + (8) x ( n ) F ( t ) = t ( n ) i + − t for t ( n ) i t t ( n ) i + . (9)These equations provide time-resolved quantities needed to define time-resolveddissimilarity profile from discrete values the same way as Eqs. 1 and 2 providethem for A-ISI-distance. The weighted spike time difference for a spike train is thencalculated as an interpolation from one difference to the next by S n ( t ) = ∆ t ( n ) i ( t ) x ( n ) F ( t ) + ∆ t ( n ) i + ( t ) x ( n ) P ( t ) x ( n ) ISI ( t ) , t ( n ) i t t ( n ) i + . (10)This continuous function is analogous to term x ( n ) ISI for the ISI-distance, except thatit is piecewise linear instead of piecewise constant. The pairwise A-SPIKE-distanceprofile is calculated by temporally averaging the weighted spike time differences,normalizing to the local firing rate average and, finally, weighting each profile bythe instantaneous firing rates of the two spike trains: S Am , n ( t ) = S n x m ISI ( t ) + S m x n ISI ( t ) h x n , m ISI ( t ) i max {h x n , m ISI ( t ) i , thr } , (11) easures of spike train synchrony and directionality 7 (a) S ( t ) (b) Time [a.u.] S A ( t ) (c) Fig. 3
An example spike train pair and its SPIKE-distance and A-SPIKE-distance profiles. (a)Two spike trains having four events with five spikes per event in each spike train. The sequence ofspikes in all four events is the same but the event is increasingly compressed. The only thing thatchanges is the time-scale. From a global perspective the first event consists of non-synchronousindividual spikes, while the last event consists of coincident bursts. The two events in the middleare intermediates. (b) The SPIKE-distance S ( t ) looks only at the local context and has the sameprofile shape for all events. (c) The A-SPIKE-distance considers also the global context and judgesthe first event like the SPIKE-distance as being dissimilar, but scales down the small spike timedifferences in the burst and considers the coincident burst as very similar. where h x n , m ISI ( t ) i is the mean over the two instantaneous ISIs. MRTS is defined byusing a threshold, that replaces the denominator of weighting to spike time differ-ences if the mean is smaller than the thr . This profile is analogous to the pairwiseA-ISI-profile I An , m ( t ) , but again it is piecewise linear, not piecewise constant. Un-like S n ( t ) it is not continuous, typically it exhibits instantaneous jumps at the timesof the spikes. The multivariate A-SPIKE-profile is obtained the same way as themultivariate A-ISI-profile, by averaging over all pairwise profiles: S A ( t ) = N ( N − ) / N − ∑ n = N ∑ m = n + S Am , n ( t ) . (12)The final A-SPIKE-distance is calculated as the time integral over the multivariateA-SPIKE-profile the same way as the A-ISI-distance: D AS = t e − t s Z t e t s S A ( t ) dt . (13)The effect of applying the threshold can be seen in Fig. 3. With thr = Eero Satuvuori and Irene Malvestio and Thomas Kreuz test dataset in Fig. 2c goes to zero when the spikes in all spike trains appear at theexactly same time.
A-SPIKE-synchronization [23] quantifies how many of the possible coincidences ina dataset are actually coincidences (Fig. 1c). While the A-ISI-distance and the A-SPIKE-distance are measures of dissimilarity which obtain low values for similarspike trains, A-SPIKE-synchronization measures similarity. If all the spikes are co-incident with a spike in all the other spike trains, the value will be one. In contrast,if none of the spikes are coincident, it will be zero.The original SPIKE-synchronization [16] is parameter- and time-scale-free, sinceit uses the adaptive coincidence detection which was first proposed for the measureEvent synchronization [21]. The coincidence window, i.e., the time lag below whichtwo spikes from two different spike trains, t ( n ) i and t ( m ) j , are considered to be coin-cident, is adapted to the local firing rate. Spikes are coincident only if they bothlie in each others coincidence windows. A-SPIKE-synchronization is a generalizedversion of the SPIKE-synchronization. The MRTS is used to decide if the windowis determined locally or if the global context should be taken into account.As a first step, we define the ISI before and after the spike as x ( n ) iP = lim t → t i − x ( n ) ISI ( t ) (14) x ( n ) iF = lim t → t i + x ( n ) ISI ( t ) . (15)The coincidence window for spike i of spike train n is defined by determining a min-imum coincidence window size for a spike as half of the ISIs adjacent to the spikeand allowing asymmetric coincidence windows based on MRTS. This is done byusing thr instead of the minimum, if it is smaller. Since the threshold value is basedon ISIs and the coincidence window spans both sides of the spike, only half of thethreshold spans each side. For the A-ISI- and the A-SPIKE-distance the changesinduced by the threshold appear gradually, but for A-SPIKE-synchronization it is asudden change from a non-coincidence to coincidence for a spike. Therefore, dueto the binary nature of A-SPIKE-synchronization, the threshold is additionally di-vided by two. The coincidence window is not allowed to overlap with a coincidencewindow of another spike and is thus limited to half the ISI even if the threshold islarger. The base of the window is defined by the two adjacent ISIs: τ ( n ) i =
12 min { x ( n ) iP , x ( n ) iF } . (16)The coincidence window of a spike is then defined in an asymmetric form by usingthe coincidence window part before and after the spike as easures of spike train synchrony and directionality 9 τ ( n ) iP = min { max ( thr , τ ( n ) i ) , x ( n ) iP } (17) τ ( n ) iF = min { max ( thr , τ ( n ) i ) , x ( n ) iF } . (18)The combined coincidence window for spikes i and j is then defined as τ ( n , m ) i j = ( min { τ ( n ) iF , τ ( m ) jP } if t i t j min { τ ( n ) iP , τ ( m ) jF } otherwise . (19)The coincidence criterion can be quantified by means of a coincidence indicator C ( n , m ) i = ( j {| t ( n ) i − t ( m ) j |} < τ ( n , m ) i j . (20)This definition ensures that each spike can only be coincident with at most one spikein the other spike train. The coincidence criterion assigns either a one or a zero toeach spike depending on whether it is part of a coincidence or not. For each spikeof every spike train, a normalized coincidence counter C ( n ) i = N − ∑ m = n C ( n , m ) i (21)is obtained by averaging over all N − i in spike train n .This way we have defined a coincidence indicator for each individual spike in thespike trains. In order to obtain one combined similarity profile, we pool the spikesof the spike trains as well as their coincidence indicators by introducing one overallspike index k and defining M = N ∑ n = M n . (22)This yields one unified set of coincidence indicators C k : { C k } = [ n { C ( n ) i } . (23)From this discrete set of coincidence indicators C k the A-SPIKE-synchronizationprofile C A ( t k ) is obtained by C A ( t k ) = C k .Finally, A-SPIKE-synchronization is defined as the average value of this profile S AC = M M ∑ k = C A ( t k ) . (24) S p i k e t r a i n S C (a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a)(a) S p i k e t r a i n S CA (b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b)(b) Time S p i k e t r a i n Diff (c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)(c)
Fig. 4
SPIKE-synchronization, A-SPIKE-synchronization and their difference illustrated usingfive spike trains with four simple events. Without the correction (A) in case of spike doublets (firstand second event) or bursts (third event) the small interspike intervals result in an unreasonablyhigh demand for spike timing accuracy. With the adaptive correction (B) for all these cases thelikelihood increases that at least one of the spikes is part of a coincidence. On the other hand, ifthere are no doublets or bursts (last event), nothing changes (best seen in C). Note that for bettervisibility the colour scales differ, we use grey-black in A and B and white-black in C.
It is important to note that since A-SPIKE-synchronization is a measure of similar-ity, reducing differences below threshold adds coincidences and thus the value ob-tained increases. In Fig. 4 we illustrate how the asymmetric coincidence windows ofA-SPIKE-synchronization allow for a better coverage of burst events which makesit easier to match spikes when compared to SPIKE-synchronization ( thr =
0) [16].As can be seen in Fig. 2d, the A-SPIKE-synchronization profile is discrete andonly defined at spike times. A dotted line between the points is added as visual aid.The profile gets higher values the more coincidences are found for each spike inother spike trains.
In some cases spikes that occur less than a second apart might be considered moresimultaneous than those taking place within minutes, and in applications like me-teorological systems, weeks instead of months. Setting the minimum relevant time- easures of spike train synchrony and directionality 11 scale might not be a simple task. If no information of the system producing thespikes is available, all one can do to estimate an appropriate threshold value is tolook at the ISIs.There are two criteria that a threshold value extracted from the data has to fulfil.First of all it needs to adapt to changes in spike count so that adding more spikesgives shorter threshold. Additionally we want the threshold to adapt to changes inthe ISI-distribution when the spike count is fixed. The more pronounced bursts arefound in the data, the more likely any differences within aligned bursts are not asimportant as their placement. Thus, we want our threshold to get longer if the spikesare packed together. To do so, all the ISIs in all spike trains are pooled and thethreshold is determined from the pooled ISI distribution.One should not just take a value based on ISI-distribution that counts the inter-spike intervals, as the mean does, but weight them by their length, which is equiv-alent to taking the average of the second moments of ISIs. Doing this reduces theimportance of very short ISIs even if they are statistically much more common. Inorder to obtain a value with the right dimension, the square root of the average mustbe taken: thr = q h ( L ISI ) i = s ∑ Nn = a n L n ISI2 ∑ Nn = a n . (25)Here we denoted a single ISI length in the pooled distribution as L n ISI and the numberof ISI with length L n ISI as a n . It is important to note, however, that this is only anestimate based on different time-scales found in the data. The selected MRTS is notan indicator of a time-scale of the system that produced the spikes.As an example of how the threshold works we apply the threshold to Gamma Γ ( k , x ) distribution. Since the kurtosis of the distribution is proportional to 1 / k , forsmall k the distribution contains large number of small ISI and few long ones.Thisis the property the threshold is tracking. The mean of a gamma distribution is k / x and the second moment ( k + ) k / x . thus the ratio of the threshold and the mean ISIis thr / h L ISI i = p x ( k + ) / k . From the formula we can see that for small k, wherethe distribution is more skewed, the ratio between the mean ISI and the thresholdincreases. This means that mainly the rare and large inter-burst ISIs are taken intoaccount.The threshold value determines the outcome of the adaptive methods. However,the threshold is not a hard set limiter neglecting everything below the threshold, butrather the point from which on differences are considered in a global instead of alocal context. Often a set of spike trains exhibits well-defined patterns of spatio-temporal prop-agation where some prominent feature first appears at a specific location and then spreads to other areas until potentially becoming a global event. If a set of spiketrains exhibits perfectly consistent repetitions of the same global propagation pat-tern, this can be called a synfire pattern . For any spike train set exhibiting propa-gation patterns the questions arises naturally whether these patterns show any con-sistency, i.e., to what extent do the spike trains resemble a synfire pattern, are therespike trains that consistently lead global events and are there other spike trains thatinvariably follow these leaders?
Spike trains (b)
Sorted spike trainsTime [a.u.] (c) (a)
Perfect synfire pattern
Fig. 5
Motivation for SPIKE-order and Spike Train Order. (a) Perfect Synfire pattern. (b) Unsortedset of spike trains. (c) The same spike trains as in (b) but now sorted from leader to follower.
In the second part of this chapter we describe a framework consisting of twodirectional measures (
SPIKE-Order and
Spike Train Order ) that allows to definea value termed
Synfire Indicator which quantifies the consistency of the leader-follower relationships [17]. This Synfire Indicator attains its maximal value of 1for a perfect synfire pattern in which all neurons fire repeatedly in a consistent orderfrom leader to follower (Fig. 5a).The same framework also allows to sort multiple spike trains from leader to fol-lower, as illustrated in Figs. 5b and 5c. This is meant purely in the sense of temporalsequence. Whereas Fig. 5b shows an artificially created but rather realistic spiketrain set, in Fig. 5c the same spike trains have been sorted to become as close aspossible to a synfire pattern. Now the spike trains that tend to fire first are on topwhereas spike trains with predominantly trailing spikes are at the bottom.Analyzing leader-follower relationships in a spike train set requires a criterionthat determines which spikes should be compared against each other. What is needed easures of spike train synchrony and directionality 13 is a match maker, a method which pairs spikes in such a way that each spike ismatched with at most one spike in each of the other spike trains. This match makeralready exists. It is the adaptive coincidence detection first used as the fundamentalingredient for the bivariate measure
SPIKE-synchronization [16] (see Section 2.3).
The symmetric measure SPIKE-Synchronization (introduced in Section 2.3) assignsto each spike of a given spike train pair a bivariate coincidence indicator. These co-incidence indicators C ( n , m ) i (Eq. 20), which are either 0 or 1, are then averaged overspike train pairs and converted into one overall profile C ( t k ) normalized between0 and 1. In exactly the same manner SPIKE-Order and Spike Train Order assignbivariate order indicators to spikes. Also these two order indicators, the asymmet-ric D ( n , m ) i and the symmetric E ( n , m ) i , which both can take the values −
1, 0, or + D ( t k ) and E ( t k ) which are normalized between − D ( t k ) dis-tinguishes 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 trainsfrom leaders to followers.First of all, the symmetric coincidence indicator C ( n , m ) i of SPIKE-Synchronization(Eq. 20) is replaced by the asymmetric SPIKE-Order indicator D ( n , m ) i = C ( n , m ) i · sign ( t ( m ) j − t ( n ) i ) , (26)where the index j is defined from the minimum in Eq. 20 with the threshold-valuein Eqs. 17 and 18 set to thr = D ( m , n ) j is obtained in an asymmetric manner as D ( m , n ) j = C ( m , n ) j · sign ( t ( n ) i − t ( m ) j ) = − D ( n , m ) i . (27)Therefore, this indicator assigns to each spike either a 1 or a − C ( n , m ) i = t ( m ) j = t ( n ) i ).The multivariate profile D ( t k ) obtained analogously to Eq. 23 is normalized be-tween 1 and − + −
1) coincident spikes in all other spike trains. It can be 0 either if aspike is not part of any coincidences or if it leads exactly as many spikes from otherspike trains in coincidences as it follows. From the definition in Eqs. 26 and 27 itfollows immediately that C k is an upper bound for the absolute value | D k | . Spike trains (a) −1010 100 200 300 400 500 600 700 800 900 1000−101 F u =−0.181(p >> 0.05) C=0.981
Time profile E (b)(c)
Pairwise matrix D
Sorted spike trainsTime (e) −1010 100 200 300 400 500 600 700 800 900 1000−101 F s =0.211 C=0.981
Time profile E for sorted spike trains (d)
Sorted pairwise matrix D F s z = −1.05 ; p >> 0.05F (f) Fig. 6
Illustration of SPIKE-Order for an artificially created example dataset consisting of 6 spiketrains which emit spikes in nine reliable events. For the first two events spikes fire in order, for thenext three events the order is random whereas for the last four events the order is inverted. In the lastevent there is one spike missing. (a) Unsorted spike trains with the spikes color-coded accordingto the value of the SPIKE-Order D ( t k ) . (b) Spike Train Order profile E ( t k ) . Events with differentfiring order can clearly be distinguished. The SPIKE-Synchronization profile C ( t k ) and its mirrorprofile (dashed black lines) act as envelope. The Synfire Indicator F u for the unsorted spike trains isslightly negative reflecting the dominance of the inversely ordered events. (c) Pairwise cumulativeSPIKE-Order matrix D before (left) and after (right) sorting. The optimal order maximizes theupper triangular matrix D ( n < m ) , marked in black. The thick black arrow in between the two matricesindicates the sorting process. (d) Spike Train Order profile E ( t k ) and its average values, the SynfireIndicator F s for the sorted spike trains. (e) Sorted spike trains. (f) Statistical significance: Resultsof the surrogate analysis. Thick lines denote mean and standard deviation for 19 surrogates. Sincethe value for the original dataset (black) is not maximum, the optimally sorted spike trains do notexhibit a statistically significant synfire pattern.easures of spike train synchrony and directionality 15 In Fig. 6 we show the application of the SPIKE-Order framework to an exampledataset. While the SPIKE-Order profile can be very useful for color-coding andvisualizing local spike leaders and followers (Fig. 6a), it is not useful as an overallindicator of Spike Train Order. The profile is invariant under exchange of spiketrains, i.e. it looks the same for all events no matter what the order of the firing is (inour example only the last event looks slightly different since one spike is missing).Moreover, summing over all profile values, which is equivalent to summing overall coincidences, necessarily leads to an average value of 0, since for every leadingspike ( +
1) there has to be a following spike ( − E ( n , m ) i = C ( n , m ) i · ( sign ( t ( m ) j − t ( n ) i ) if n < m sign ( t ( n ) i − t ( m ) j ) if n > m (28)and E ( m , n ) j = E ( n , m ) i . (29)This symmetric indicator assigns to both spikes a + − E ( t k ) , again obtained similarly to Eq. 23, is also normal-ized between 1 and − C k is anupper 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 itdoes make sense to define an average value, which we term the Synfire Indicator: F = M M ∑ k = E ( t k ) . (30) The interpretation is very intuitive. The Synfire Indicator F quantifies to whatdegree the spike trains in their current order resemble a perfect synfire pattern. Itis normalized between 1 and − −
1) if the spike trains intheir current order form a perfect (inverse) synfire pattern, meaning that all spikesare coincident with spikes in all other spike trains and that the order from leading(following) to following (leading) spike train remains consistent for all of theseevents. It is 0 either if the spike trains do not contain any coincidences at all or ifamong all spike trains there is a complete symmetry between leading and followingspikes.The Spike Train Order profile E ( t k ) for our example is shown in Fig. 6c. In thiscase the order of spikes within an event clearly matters. The Synfire Indicator F isslightly negative indicating that the current order of the spike trains is actually closerto an inverse synfire pattern.Given a set of spike trains we now would like to sort the spike trains from leaderto follower such that the set comes as close as possible to a synfire pattern. To doso we have to maximize the overall number of correctly ordered coincidences andthis is equivalent to maximizing the Synfire Indicator F . However, it would be verydifficult 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 ofthe spike trains. The most intuitive way is to use the anti-symmetric cumulativeSPIKE-Order matrix D ( n , m ) = ∑ i D ( n , m ) i (31)which sums up orders of coincidences from the respective pair of spike trains onlyand quantifies how much spike train n is leading spike train m (Fig. 6c).Hence if D ( n , m ) > n is leading m , while D ( n , m ) < m isleading n . If the current Spike Train Order is consistent with the synfire property,we thus expect that D ( n , m ) > n < m and D ( n , m ) < n > m . Therefore, weconstruct the overall SPIKE-Order as D ( n < m ) = ∑ n < m D ( n , m ) , (32)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 asecond more practical definition of the Synfire Indicator: F = D ( n < m ) ( N − ) M . (33)The value is identical to the one of Eq. 30, only the temporal and the spatial summa-tion of coincidences (i.e., over the profile and over spike train pairs) are performedin the opposite order.Having such a quantification depending on the order of spike trains, we can in-troduce a new ordering in terms of the spike train index permutation ϕ ( n ) . Theoverall Synfire Indicator for this permutation is then denoted as F ϕ . Accordingly, easures of spike train synchrony and directionality 17 for the initial ( u nsorted) order of spike trains ϕ u the Synfire Indicator is denoted as F u = F ϕ u .The aim of the analysis is now to find the optimal ( s orted) order ϕ s as the oneresulting in the maximal overall Synfire Indicator F s = F ϕ s : ϕ s : F ϕ s = max ϕ { F ϕ } = F s . (34)This Synfire Indicator for the sorted spike trains quantifies how close spike trainscan be sorted to resemble a synfire pattern, i.e., to what extent coinciding spike pairswith correct order prevail over coinciding spike pairs with incorrect order. Unlike theSynfire 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 couldsimply be reversed in order to obtain the same positive value). For a perfect synfirepattern we obtain F s =
1, while sufficiently long Poisson spike trains without anysynfire structure yield F s ≈ N spike trains there are N !permutations ϕ , so for large numbers of spike trains finding the optimal Spike TrainOrder ϕ s is a non-trivial problem and brute-force methods such as calculating the F ϕ -value for all possible permutations are not feasible. Instead, we search for theoptimal order using simulated annealing [7], a probabilistic technique which ap-proximates the global optimum of a given function in a large search space. In ourcase 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 usingthe fundamental move of exchanging two neighboring spike trains within the cur-rent permutation. All moves with positive ∆ F are accepted while the likelihood ofaccepting moves with negative ∆ F is decreased along the way according to a stan-dard slow cooling scheme. The procedure is repeated iteratively until the order ofthe spike trains no longer changes or until a predefined end temperature is reached.The sorting of the spike trains maximizes the Synfire Indicator as reflected byboth the normalized sum of the upper right half of the pairwise cumulative SPIKE-Order matrix (Eq.33, Fig. 6c) and the average value of the Spike Train Order profile E ( t k ) (Eq.30, Fig. 6d). Finally, the sorted spike trains in Fig. 6e are now orderedsuch that the first spike trains have predominantly high values (red) and the lastspike 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 theprofiles D and E . The pairwise information in the SPIKE-order matrix reflects theleader-follower relationship between two spike trains at a time. The Synfire Indica-tor F characterizes the closeness of the dataset as a whole to a synfire pattern, bothfor the unsorted ( F u ) and for the sorted ( F s ) spike trains. Finally, the sorted order ofthe spike trains is a very important result in itself since it identifies the leading andthe following spike trains. As a last step in the analysis we evaluate the statistical significance of the opti-mized Synfire Indicator F s using a set of carefully constructed spike train surrogates.The idea behind the surrogate test is to estimate the likelihood that the consistentSPIKE-Order pattern yielding a certain Synfire Indicator could have been obtainedby chance. To this aim, for each surrogate we maintain the coincidence structureof the spike trains by keeping the SPIKE-Synchronization values of each individualspike constant but randomly swap the spike order in a sufficient number of coin-cidences. We set the number of swaps equal to the number of coincident spikes inthe dataset since this way all possible spike order patterns can be reached. Only forthe first surrogate we swap twice as many coincidences in order to account for tran-sients. After each swap we take extra care that all other spike orders that are affectedby the swap are updated as well. For example, if a swap changes the order betweenthe first and the third spike in an ordered sequence of three spikes, we also swapboth the order between the first and the second and the order between the secondand the third spike.For each spike train surrogate we repeat exactly the same optimization procedurein the spike train permutation space that is done for the original dataset. The origi-nal Synfire Indicator is deemed significant if it is higher than the Synfire Indicatorobtained 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 ∗ = / ( s + ) = . s is acompromise that takes into account the computational cost. As a second indicatorwe state the z-score, e.g., the deviation of the original value x from the mean µ ofthe surrogates in units of their standard deviation σ : z = x − µσ . (35)Results of the significance analysis for our standard example are shown in thehistogram in Fig. 6f. In this case the absolute value of the z-score is smaller thanone and the p-value is larger than p ∗ and the result is thus judged as statisticallynon-significant. In the first part of this chapter we describe three parameter-free and time resolvedmeasures of spike train synchrony in their recently developed adaptive extensions,A-ISI-distance, A-SPIKE-distance and A-SPIKE-synchronization [23]. All of thesemeasures are symmetric and so their multivariate versions are invariant to changesin the order of spike trains. Since information about directionality is very relevant, easures of spike train synchrony and directionality 19 in the second part of this chapter we show an algorithm which allows to sort mul-tiple spike trains from leader to follower. This algorithm is built on two indicators,SPIKE-Order and Spike Train Order, that define the Synfire Indicator value, whichquantifies the consistency of the temporal leader-follower relationships for both theoriginal and the optimized sorting.Symmetric measures of spike train distances (Section 2) have been applied inmany different contexts, not only in the field of neuroscience [10, 18, 27]. For exam-ple, they have been used in robotics [9] and prosthesis control [8]. The ISI-distancehas been applied in a method for the detection of directionl coupling between pointprocesses and point processes and flows [1], as an adaptation of the nonlinear tech-nique for directional coupling detection of continuous signals [5].Questions about leader-follower dynamics (Section 3) have been specifically in-vestigated in neuroscience [20], but also in fields as wide-ranging as, e.g., clima-tology [3], social communication [26], and human-robot interaction [22]. SPIKE-Order has already been applied to analyze the consistency of propagation patternsin two real datasets from neuroscience (Giant Depolarized Potentials in mice slices)and climatology (El Ni˜no sea surface temperature recordings) [17].Finally, we would like to mention that the similarity measures A-ISI-distance,A-SPIKE-distance and A-SPIKE-synchronization, as well as SPIKE-Order, are im-plemented in three publicly available software packages, the Matlab-based graphicaluser interface SPIKY [16], cSPIKE (Matlab command line with MEX-files), andthe open-source Python library PySpike [19]. Acknowledgements
We acknowledge funding from the European Union’s Horizon 2020 researchand innovation program under the Marie Skodowska-Curie Grant Agreement No.
References
1. R. G. Andrzejak and T. Kreuz. Characterizing unidirectional couplings between point pro-cesses and flows.
Europhysics Letters , 96:50012, 2011.2. D. L. Applegate, R. E. Bixby, V. Chvatal, and W. J. Cook.
The traveling salesman problem: acomputational study . Princeton University Press, 2011.3. N. Boers, B. Bookhagen, H.M.J. Barbosa, N. Marwan, J. Kurths, and J.A. Marengo. Predictionof extreme floods in the eastern central andes based on a complex networks approach.
NatureCommunications , 5:5199, 2014.4. M. R. Bower, M. Stead, F. B. Meyer, W. R. Marsh, and G. A. Worrell. Spatiotemporal neuronalcorrelates of seizure generation in focal epilepsy.
Epilepsia , 53:807, 2012. http://mariomulansky.github.io/PySpike0 Eero Satuvuori and Irene Malvestio and Thomas Kreuz5. D. Chicharro and R. G. Andrzejak. Reliable detection of directional couplings using rankstatistics. Physical Review E , 80(2):026217, 2009.6. D. Chicharro, T. Kreuz, and R. G. Andrzejak. What can spike train distances tell us about theneural code?
J Neurosci Methods , 199:146–165, 2011.7. K. A. Dowsland and J. M. Thompson. Simulated annealing. In
Handbook of Natural Com-puting , pages 1623–1655. Springer, 2012.8. S. Dura-Bernal, K. Li, S. A. Neymotin, J. T. Francis, J. C. Principe, and W. W. Lytton. Restor-ing behavior via inverse neurocontroller in a lesioned cortical spiking model driving a virtualarm.
Frontiers in Neuroscience , 10, 2016.9. A. Espinal, H. Rostro-Gonzalez, M. Carpio, E. I. Guerra-Hernandez, M. Ornelas-Rodriguez,H. J. Puga-Soberanes, M. A. Sotelo-Figuero, and P. Melin. Quadrupedal robot locomotion: abiologically inspired approach and its hardware implementation.
Computational Intelligenceand Neuroscience , page 5615618, 2016.10. R. Jolivet, R. Kobayashi, A. Rauch, R. Naud, S. Shinomoto, and W. Gerstner. A benchmarktest for a quantitative assessment of simple neuron models.
J Neurosci Methods , 169:417,2008.11. T. Kreuz. Synchronization measures. In R. Quian Quiroga and S. Panzeri, editors,
Principlesof neural coding , page p. 97. CRC Taylor and Francis, Boca Raton, FL, USA, 2013.12. T. Kreuz, D. Chicharro, M. Greschner, and R. G. Andrzejak. Time-resolved and time-scaleadaptive measures of spike train synchrony.
J Neurosci Methods , 195:92, 2011.13. T. Kreuz, D. Chicharro, C. Houghton, R. G. Andrzejak, and F. Mormann. Monitoring spiketrain synchrony.
J Neurophysiology , 109:1457, 2013.14. T. Kreuz, J. S. Haas, A. Morelli, H. D. I. Abarbanel, and A. Politi. Measuring spike trainsynchrony.
J Neurosci Methods , 165:151, 2007.15. T. Kreuz, F. Mormann, R. G. Andrzejak, A. Kraskov, K. Lehnertz, and P. Grassberger. Measur-ing synchronization in coupled model systems: A comparison of different approaches.
PhysD , 225:29, 2007.16. T. Kreuz, M. Mulansky, and N. Bozanic. SPIKY: A graphical user interface for monitoringspike train synchrony.
J Neurophysiol , 113:3432, 2015.17. Thomas Kreuz, Eero Satuvuori, Martin Pofahl, and Mario Mulansky. Leaders and followers:Quantifying consistency in spatio-temporal propagation patterns.
New Journal of Physics ,19:043028, 2017.18. Z. Mainen and T. J. Sejnowski. Reliability of spike timing in neocortical neurons.
Science ,268:1503, 1995.19. M. Mulansky and T. Kreuz. Pyspike - A python library for analyzing spike train synchrony.
Software X , 5:183–189, 2016.20. E. Pereda, R. Quian Quiroga, and J. Bhattacharya. Nonlinear multivariate analysis of neuro-physiological signals.
Progress in Neurobiology , 77:1, 2005.21. R. Quian Quiroga, T. Kreuz, and P. Grassberger. Event synchronization: A simple and fastmethod to measure synchronicity and time delay patterns.
Phys. Rev. E , 66:041904, 2002.22. F. Rahbar, S. Anzalone, G. Varni, E. Zibetti, S. Ivaldi, and M. Chetouani. Predicting extraver-sion from non-verbal features during a face-to-face human-robot interaction.
InternationalConference on Social Robotics , page 10, 2015.23. E. Satuvuori, M. Mulansky, N. Bozanic, K. Lenk, and T. Kreuz. Generalizing definitions ofISI-distance, SPIKE-distance and SPIKE-synchronization.
In preparation , 2017.24. W. Truccolo, J. A. Donoghue, L. R. Hochberg, E. N. Eskandar, J. R. Madsen, W. S. Anderson,E. N. Brown, E. Halgren, and S. S. Cash. Single-neuron dynamics in human focal epilepsy.
Nature Neurosci , 14:635, 2011.25. M. C. W. van Rossum. A novel spike distance.
Neural Computation , 13:751, 2001.26. G. Varni, G. Volpe, and A. Camurri. A system for real-time multimodal analysis of nonverbalaffective social interaction in user-centric media.
IEEE Transactions on multimedia , 12:576,2010.27. J. D. Victor. Spike train metrics.
Current Opinion in Neurobiology , 15:585, 2005.28. J. D. Victor and K. P. Purpura. Nature and precision of temporal coding in visual cortex: Ametric-space analysis.