### Source Code - SPIKE- and ISI-Distance

__SPIKE- and ISI-Distance__ between two or more
spike trains

For a

**detailed description**of the methods please refer to:

Kreuz T, Chicharro D, Houghton C, Andrzejak RG,
Mormann F:

Monitoring spike train synchrony.

*J Neurophysiol* **109**, 1457
(2013) [PDF].

and

Kreuz T:

SPIKE-distance.

*Scholarpedia* **7**(12), 30652 (2012).

**overview**please refer to:

Mulansky M, Bozanic N, Sburlea A, Kreuz T:

A guide to time-resolved and parameter-free measures of spike train synchrony.

IEEE Proceeding on Event-based Control, Communication, and Signal Processing (EBCCSP), 1-8 and arXiv [PDF] (2015).

**Powerpoint presentation**

Kreuz T:

Monitoring spike train synchrony.

[PDF].

Further information on the SPIKE-distance can also be found on the SPIKY Facebook page.

The best
demonstration of the SPIKE-distance can be seen in the following movie.
This and many more movies can now also be seen on the SPIKY Youtube Channel.** **

Movie: How the SPIKE-distance traces spike train clustering

Here we explain how this
time-resolved and parameter-free measure traces spike train clustering
with both its normal and its realtime version.

Top: Rasterplot with artificially generated spike trains whose clustering changes every 500 ms.

Middle: Dissimilarity matrices for the regular (1st column) and the
real-time (2nd column) SPIKE-distance. The spike train group matrices
(3rd and 4th column) are obtained by averaging over the respective
submatrices of the original dissimilarity matrices.

Bottom: Corresponding dendrograms.

In the first part of the movie the instantaneous dissimilarities and
the clustering dendrograms of both measures are updated as the green
time bar moves from left to right. The second part comprises several
examples of selective temporal averaging of individual or combined
intervals as well as triggered averaging.

A second example:

**Movie: Demonstration of externally triggered averaging.**

Here we explain how the instantaneous
nature of the SPIKE-distance can be very useful when this method is
combined with triggered averaging.

In this simulated setup a population of 50 neurons is recorded upon presentation of a time-varying non-periodic stimulus, in this case a chirp function. It is assumed that half of the neurons are sensitive to negative amplitudes and accordingly exhibit higher (lower) reliability for local minima (maxima) of the chirp function.

As the amplitude of the chirp function varies, so does the spike train synchrony of half of the neurons, but due to the non-periodicity of the stimulus this is quite difficult to detect. Detection is facilitated by externally triggered averaging where the triggering is performed on common stimulus features, here the amplitude of the chirp function. As can be seen in the dissimilarity matrices of the SPIKE-distance and even better in the dendrogram, a decrease of this stimulus amplitude leads to the emergence of a spike train cluster consisting of the modulated spike trains which indicates their increase in reliability.

The definiton of the SPIKE-distance builds on these earlier papers:

Kreuz T, Chicharro D, Greschner M, Andrzejak RG:

Time-resolved and time-scale adaptive measures of spike train synchrony.

*J Neurosci Methods* **195**, 92 (2011) [PDF].

Kreuz T, Chicharro D, Andrzejak RG, Haas JS, Abarbanel HDI:

Measuring multiple spike train synchrony.

*J Neurosci Methods* **183**, 287 (2009) [PDF].

Kreuz T, Haas JS, Morelli A, Abarbanel HDI, Politi A:

Measuring spike train synchrony.

*J Neurosci Methods* **165**, 151 (2007) [PDF].

The Matlab code is now complemented by PySpike, an open source Python library written by Mario Mulansky and available on github. Its core functionality is the implementation of the bivariate ISI- and SPIKE-distance. Additionally, it provides functions to compute multivariate ISI- and SPIKE-distances, as well as averaging and general spike train processing. All computation intensive parts are implemented in C (via cython) to reach a competitive performance (factor 100-200 over plain Python). There is also documentation, currently consisting of a brief introduction and the API reference.

See also:

Python-Implementation of the pairwise SPIKE-distance (written by Jeremy Fix)

Python-Implementation of the pairwise ISI-distance (maintained by Michael Chary)

===================================

**Matlab code to calculate and visualize the SPIKE- and the ISI-Distance between two or more
spike trains:**

**
**Kreuz_et_al_SPIKE-distance.zip
**
**

(**Last updated version** 1.2: April 29,
2013; **Copyright**: Thomas Kreuz, Nebojsa Bozanic)

**Important note:** This code will not be updated anymore. In principle it is replaced by the program 'SPIKY_loop' which is included in the zip-file of the new **graphical user interface** **SPIKY** which makes use of the same MEX-files. The Beta-version of these source codes can now be found
here.

This version uses **MEX**-
instead of m-files for the most time-consuming calculations. This reduces
the computational cost considerably (on my notebook the code running all the
examples below was faster by a factor of 85). It also uses a new **output-structure**
**'results' **which allows easy access not only to the overall
dissimilarity value but also to the measure profiles and the pairwise
distance matrices for all the selected measures.

**Note of warning (see also 'Notice' at the end):**

This program calculates the SPIKE- and the ISI-Distance between two or more spike trains. It also allows the time-resolved visualization via measure profiles and the spatial visualization via dissimilarity matrices (plots and figures). This, however, makes it a rather complex program so it might happen initially that there are some bugs or that the program is not as error-tolerant as it should be. Whenever you encounter any problem please provide feedback and we will try to resolve the problem and improve the source codes accordingly.

The name of the zip-file
will always include the date of release which will allow users to ensure
that they have the most recent version. As the development of the measure evolved, more and more
applications came to mind and were integrated, and this might not be the end
of it. Thus it might be worthwhile to check for updates now and then, in
particular in these early alpha development stages. If you want to be
informed about possible future releases, please send an email to Thomas Kreuz at "**thomas.kreuz (at) cnr.it**". You will then be included in a
mailing list.

**Description:**

The zip-file contains many m- and eight MEX-files which
all should be stored in the same
directory. Before the first run of the program the MEX-files should be
compiled. To do so please run from this directory the program **'Compile_MEX.m'**.

**Distances_Main_Demo**: Demo program, used to create all the figures below

**f_get_data****_demo**: Main
function (sets spike-matrix as well as the parameters describing these
data), called by 'Distances_Main_Demo'

**Distances_Main**: Main program; loads or creates data, sets
parameters, and calls main function 'f_distances'

**f_get_data**: Main function (sets spike-matrix as well as the
parameters describing these data), called by 'Distances_Main'

**f_distances_MEX**: Main function (calculates distances and does the
plotting), called by 'Distances_Main' or 'Distances_Main_Demo',
uses MEX-functions for computationally expensive loops

This function calls the following functions (included):

**f_measure_profiles**: Function that creates individual measure
profiles, called by 'f_distances_MEX'

**f_moving_average (& variants)**: Functions to generate the
(weighted) moving average of a measure profile, depending on the measure
profile the non-causal or the causal variant is used

**f_compute_gauss_smooth**: Function that creates a smoothed version
of the PSTH

(Copyright: Jude Mitchell, http://www.snl.salk.edu/~jude/sfn2008/index.html)

as well as a few more functions (some of them from Matlab toolboxes)

To get started, try either **Distances_Main_Demo **or** Distances_Main**.**
**However, before the first start of the program a few variables (see also
below) defining the work environment should be set:

**f_para.imagespath**: Path where images (ps) will be stored

**f_para.moviespath**: Path where movies (avi) will be stored

**f_para.matpath**: Path where Matlab files
(mat) will be stored

Main function **
f_distances_MEX:**

**Input:**

Matrix **spikes** with two or more spike trains (if trains have different numbers of spikes, fill with zeros)

Structure **d_para **that describe the data (see below)

Structure **
f_para **that determines the appearance of the figures and the movie

__
Output:__

Structure** 'results' **with these
fields:

all of these variables have one entry for
each measure selected via the variable '**f_para.subplot_posi**' (see
below):

**measures**: Name of selected measures
(helps to identify the order within all other variables)

**dissimilarity_profiles**: cell with
overall dissimilarity profiles obtained by averaging over spike train pairs
(either piecewise constant (size 'num_isi') or sampled (length 'len'),
see additional figure 1 below for an example)

**distance_matrices**: pairwise distance
matrices, obtained by averaging over time

**overall_dissimilarities**: overall
dissimilarities, obtained by averaging over both spike trains and time

In case a sampled measure was selected:

**time**: sampling values of length 'len'
(defines x-axis for sampled measures)

In case a piecewise constant measure was selected:

**isi**: length of successive isi within
the pooled spike train, length 'num_isi' (defines x-axis for piecewise
constant measures)

For the piecewise constant profiles (here
for the first) the function **f_pico** can be used to obtain the average
value as well as x- and y-vectors for plotting:

[overall_dissimilarity,plot_x_values,plot_y_values] = f_pico(results.isi,results.dissimilarity_profiles{1},d_para.dts,d_para.tmin);

**Example call:**

After all other parameters have been set to standard values, i.e.

d_para=d_para_default;

f_para=f_para_default;

d_para.tmin=0; d_para.tmax=101; d_para.dts=0.001;

num_trains=4; num_spikes=100;

spikes=zeros(num_trains,num_spikes);

spikes(1,1:num_spikes)=1:num_spikes;

for rc=2:num_trains

spikes(rc,1:num_spikes-rc)=spikes(1,1:num_spikes-rc)+(rand(1,num_spikes-rc)-0.5);

end

results=f_distances_MEX(spikes,d_para,f_para);

**Parameters:**

The program uses three sets of parameters : Each of these are stored in a Matlab-struct.

1. **d_para**: Parameters that describe the data (set
within the function **f_get_data** and then used in
the call of the main function **f_distances_MEX **from the main program)

2. **f_para**: Parameters that determine the
appearance of the figures and the movie (set in the main program and then used in the call of the main
function **f_distances_MEX **from the main program)

3. **s_para**: Parameters that determine the appearance of
the individual subplots in the figure (set in the function **f_distances_MEX **and then used in** **the call of the
function **f_measure_profiles**)

Here is a description of the individual structures:

1. Structure '__ d_para__': parameters that
describe the data

**tmin**: Beginning of recording

**tmax**: End of recording

**dts**: Sampling interval, precision of spike times [!!! Please take care that this value is not larger than the actual sampling size, otherwise two spikes can occur at the same time instant and this can lead to problems in the algorithm !!!]

**all_train_group_names**: Names of spike train groups

**all_train_group_sizes**: Sizes of respective spike train groups

**select_train_mode**: Selection of spike trains (1-all,2-selected groups,3-selected trains)

**select_train_groups**: Selected spike train groups (if 'select_train_mode==2')

**select_trains**: Selected spike trains (if 'select_train_mode==3')

**select_averages**: One or more continuous intervals for selective temporal averaging

**trigger_averages**: Non-continuous time-instants for triggered temporal averaging, external (e.g. certain stimulus properties) but also internal (e.g. certain event times)

**markers**: Relevant time instants

**markers2**: Even more relevant time instants

**separators**: Relevant seperations between groups of spike trains

**separators2**: Even more relevant seperations between groups of spike trains

**interval_separators**: Edges of subsections

**interval_strings**: Captions for subsections

**comment_string**: Comment on the data, will be used in
file and figure names

2. Structure '** f_para**': parameters that
determine the appearance of the figures (and the movie)

**imagespath**: Path where images (postscript) will be stored

**moviespath**: Path where movies (avi) will be stored

**matpath**: Path where Matlab file (mat) will be stored

**filename**: Name under which images, movies and Matlab files will be stored

**title_string**: Appears in the figure titles

**saving**: Saving to postscript file? (0-no,1-yes)

**print_mode**: Saving to postscript file? (0-no,1-yes)

**publication**: Omits otherwise helpful information, such as additional comments (0-no,1-yes)

**comment_string**: Additional comment on the example, will be used in file and figure names

**num_fig**: Number of figure

**pos_fig**: Position of figure

**font_size**: Font size of labels (and headlines)

**multi_figure**: Open many figures (0-no,1-yes)

**timeunit_string**: string with time unit, used in labels

**xfact**: Conversion of time unit

**ma_mode**: Moving average mode: (0-no,1-only,2-both)

**spike_mao**: Order of moving average (pooled ISI)

**time_mao**: Order of moving average (time)

**dtm**: Sampling of measure profile, setting this to larger values allows downsampling in order to facilitate memory management

**mov_step**: Step size for movie frames

**mov_frames_per_second**: Well, frames per second

**mov_num_average_frames**: Number of mov_frames the averages are shown at the end of the movie (if this is too small they are hardly visable)

**mov_frames**: Selection of individual frames to be shown in the movie (replaces use of mov_step as step size if non-empty)

**plot_mode**: +1:vs time,+2-different measures and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie (binary addition allows all combinations)

**norm_mode**: normalization of averaged bivariate measure profiles (1-Absolute maximum value 'one',2-Overall maximum value,3-Individual maximum value)

**color_norm_mode**: normalization of pairwise color matrices (1-Absolute maximum value,2-Overall maximum value,3-Each one individually)

**block_matrices**: Allows tracing the overall synchronization among groups of spike trains (0-no,1-yes)

**dendrograms**: Cluster trees from distance matrices (0-no,1-yes)

**dendro_color_mode**: coloring of dendrograms (0-no,1-groups,2-trains)

**subplot_size**: Vector with relative size of subplots

**subplot_posi**: Vector with order of subplots

3. Structure '__ s_para__': parameters that
describe the appearance of the individual subplots in the figures (measure
time profiles)

**window_mode**: time interval selection (1-all (recording limits),2:outer spikes,3-inner spikes,4-smaller analysis window)

**colors**: Colors for 1:Stim/STs/ISI,2:Stim/STs/ISI-outside,3:Mean values,4:Mean values-outside,5:Moving averages,6:Moving averages-outside

**font_size**: Font size of labels (and headlines)

**line_width**: Line width

**line_style**: Line style (1-,2:,3-.,4--,5:none)

**nma**: Selected moving averages (depends on f_para.ma_mode)

**spike_mao**: Order of moving average (pooled ISI)

**time_mao**: Order of moving average (time)

**causal**: determines kind of moving average, set automatically for each measure (0-no,1-yes)

**itmin**: Beginning of recording

**itmax**: End of recording

**wmin**: Beginning of selected window (if window_mode==4)

**wmax**: End of selected window (if window_mode==4)

**num_subplots**: depends on f_para.subplot_posi, set automatically

**xl**: x-limits for plotting, set automatically

**yl**: y-limits for plotting, set automatically

**dtm**: equals f_para.dtm

**plot_mode**: equals f_para.plot_mode

Some of the more important parameters are here explained in more detail:

**f_para.plot_mode** - This variable selects the types of plots
(and movie)

+1: time-resolved measure profile (see Figs. 2 and 3 of the 2013 paper and below)

+2: different measures and cuts (see Figs. 5 and 6 of the 2013 paper and below)

+4: different measures (creates one figures of different cuts for each measure)

+8: different cuts (basically screenshots of the movie, see Fig. 9 in the 2013 paper and below)

+16: different cuts-Movie (basically a succession of plots, see Fig. 9 in the 2013 paper and below as well as the supplementary movie)

Binary addition allows all combinations of these plots/movies.

For the __ measure profiles__ (created by adding +1 to

**f_para.plot_mode**):

**f_para.subplot_posi **
and **f_para.subplot_size - **These variables allow you
to choose the selection and order the subplots as well as their relative
size.

Use **f_para.subplot_posi **to select the order of the subplots. Use 0 if a subplot is not needed. Make sure that the range from 1 to the
desired number of subplots is covered, no (positive) number should appear twice.

Use **f_para.subplot_size **to select the relative size of the selected subplots.

The first variable **f_para.subplot_posi** refers to the standard
order of the subplots which is defined below.

For each of the following variables there will be an entry in the output variable (in the order of the list below.

The second variable **f_para.subplot_size** then refers to the subplots selected by "f_para.subplot_posi" and gives their
relative sizes from top to bottom (its length should be equal to the maximum number in "f_para.subplot_posi").

The subplot "Stimulus" allows you to add any representation of the stimulus that triggered the spike trains (put 0 if not needed).

In **f_distances_MEX** the function **f_measure_profiles** is called
once for each selected measure profile. This is a list of the measures
provided.

Stimulus, Spike trains, Classic, Rate, ISI:

1
2
3
4 5
6 7
8
9 10

Stimulus Spikes PSTH GPSTH ISI
<ISI> Rate <Rate> Ia I^a

SPIKE:

11 12
13
14
15
16
17
18 19
20
21 22

Sa
S^a Sa_r
S_r^a Sa_f
S_f^a Sta
St^a Sta_r
St_r^a Sta_f St_f^a

From 9 to 22 these are six pairs where the first entry (Xa) refers to all pairwise measure profiles and the second entry (X^a) to the averaged bivariate measure profiles.

**Piecewise constant and sampled variants**:

The dissimilarity profile S (t) of the SPIKE-distance is piecewise linear rather than piecewise constant (as is the case for the ISI-distance). Therefore, when the localized visualization is desired, a new value has to be calculated for each sampling point and not just once per each interval in the pooled spike train.

In cases where the distance value itself is sufficient, the short computation time can be even further decreased by representing each interval by the value of its center and weighting it by its length. This is not only faster, but it actually gives the exact result, whereas the time-resolved calculation is a very good approximation only for sufficiently small sampling intervals dt (imagine the example of a rectangular function, at some point any sampled representation has to cut the right angle). The dissimilarity profile S_r (t) of the real-time SPIKE-distance is hyperbolic and not linear but here also the exact result can be obtained by piecewise integration over all intervals of the pooled spike train.

Measures **9-16** are **piecewise constant** (one value for each
interval in the pooled spike train), while measures **17-22** are **sampled** with
the sampling interval given by the variable **f_para.dtm**.

9-10: ISI-distance (per definition piecewise constant)

11-12: SPIKE-distance (piecewise constant)

13-14: Real-time SPIKE-distance (piecewise constant)

15-16: Future SPIKE-distance (piecewise constant)

17-18: SPIKE-distance (sampled with variable 'time')

19-20: Real-time SPIKE-distance (sampled with variable 'time')

21-22: Future SPIKE-distance (sampled with variable 'time')

For each selected variable a value in the variable 'mean_values'
is returned (the order is defined by **f_para.subplot_posi)**.

There exist a number of different possibilities how the first and the last inter-spike interval are defined. For longer spike trains with sufficient statistics these differences become rather irrelevant.

**s_para .window_mode** (1-all, 2-outer spikes, 3-inner spikes,
4-chosen analysis window)

*1. All* - time profiles are calculated within this interval (which should cover the spike train such that in each spike
train the first interval is from **d_para.tmin** to the first spike, the last interval is from the last spike to **
d_para.tmax)**. This corresponds to
the case where spike trains were recorded for a certain time. The variable d_para.min then denotes the onset of the recording,
d_para.tmax its end.

*2. Outer spikes* - time profiles are averaged from the very first spike till the very last spike.

*3. Inner spikes*
- time profiles are calculated on the common support, i.e., from the
latest first spike till the earliest last spike. This option does not
make much sense if the spike trains do not have overlapping support
(i.e., if one spike train starts after the other has already ended).

*4. Window* - restricts the averaging to a smaller window [**s_para.wmin
s_para.wmax**] to be defined. These values should be within the range [tmin tmax].

**s_para.wmin** and **s_para.wmax** - beginning and end of the analysis window
(for window_mode=4).

**s_para.colors** - colors used to display the stimulus/spike trains/ISI,
the mean values and the moving averages both for the whole time profile and
the selected analysis window (for window_mode>1)

**f_para.publication** (0-no,
1-yes) - If yes, the figure will have a title with all the selected distance
values and the postscript-file will be a large landscape print, otherwise no
title and a smaller portrait file.

For the __ movies__ (created by adding +16 to

**f_para.plot_mode**):

The movies will always append individual figures in this order:

1. **Individual frames** (either individually selected or sampled with
a certain step size)

2. **
Selected averages**

3. **
Triggered averages**

There are two possibilities to select **individual frames** in the
first part of the movie.

1. **f_para.mov_frames** - Selection of individual frames to be shown
in the movie (replaces use of mov_step as step size if non-empty)

2. **f_para.mov_step ** - Step size for movie frames
(will only be used if f_para.mov_frames is empty).

**f_para.mov_num_average_frames** - Number of frames
the selected and triggerad averages are shown at the end of the movie (in
case the first part of the movie contains a sequence of snapshots with high
temporal resolution and f_para.mov_frames_per_second is rather large this
number has to be rather large (e.g., of the order of
f_para.mov_frames_per_second) in order to make the averages visible)

**f_para.block_matrices**: Allows tracing the overall synchronization
among groups of spike trains (0-no,1-yes) (see Fig. 9 of the 2013 paper and
below)

**f_para.dendrograms**: Cluster trees from distance matrices
(0-no,1-yes), (see Figs. 8 and 9 of the 2013 paper and below)

**d_para.select_averages**: Continuous intervals for
selective temporal averaging (see Fig. 7 of the 2013 paper and
below)

**d_para.trigger_averages**: Non-continuous
time-instants for triggered temporal averaging, external (e.g. certain
stimulus properties) but also internal (e.g. certain event times) (see Fig.
8 of the 2013 paper and below)

To control the position and size of the matrix or
dendrogram subplots the variable **supos **is set in **f_distances_MEX**.
Standard settings for up to 8 different subplots are provided, however,
further adjustments might be necessary.

In case you want to create your own
plots the pairwise results are stored in the four-dimensional variable '**mov_mat**'
with dimensions [#(selected measures)] * [#(selected frames) + #(selected
averages) + #(selected triggers)] *[#(spike trains)[ * [#(spike trains)].
The respective variables for the second dimensions are num_frame_select,
num_select_averages, and num_trigger_averages.

**Memory management**:

For large numbers of spike trains with many spikes the interval is cut into smaller segments, and the averaging over all pairs of spike trains is performed for each segment separately. During calculations this is indicated on the screen by 'run_info = #run #(overall runs)'.

In the beginning of the function the
variable '**max_memo_init**' is set. This variable should be big enough
to hold the basic matrices for one run of the segment loop but small enough
to not run out of memory.

**===================================================**

In the following you will find an explanation of the most relevant parameter settings needed in order to reproduce some of the plots from the 2013 paper:

To reproduce these figures run 'Distance_Main_Demo' and select the figures by setting the variable 'examples' accordingly

**Fig. 2a: SPIKE-distance applied to a bivariate
example with varying spike match**

(example 21)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0
0 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.comment_string=['Fig2a']; % Additional comment on the example, will
be used in file and figure names

**Fig. 2b: SPIKE-distance applied to a multivariate
example**

(example 22)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0
0 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.comment_string=['Fig2b']; % Additional comment on the example, will
be used in file and figure names

**Fig. 3a: Real-time SPIKE-distance applied to a bivariate example with varying spike match**

(example 31)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 2 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.ma_mode=2; % Moving average mode: (0-no,1-only,2-both)

f_para.time_mao=120; % Order of moving average (time)

f_para.comment_string=['Fig3a_Realtime']; % Additional comment on the
example, will be used in file and figure names

**Fig. 3b: SPIKE-distance applied to a multivariate
example**

(example 32)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 2 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.ma_mode=2; % Moving average mode: (0-no,1-only,2-both)

f_para.time_mao=190; % Order of moving average (time)

f_para.comment_string=['Fig3b_Realtime']; % Additional comment on the
example, will be used in file and figure names

**Fig. 4: Real-time SPIKE-distance: Peaks during
reliable spiking events
are not spurious**

(example 41)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 2 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.comment_string=['Fig4']; % Additional comment on the example, will be
used in file and figure names

**Fig. 5A: Instantaneous clustering for artificially
generated spike trains**

(example 51)

d_para.interval_separators=500:500:d_para.tmax-500;
% Edges
of subsections

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 3 0 0]; % Vector with order of subplots

f_para.plot_mode=2; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.mov_frames=[250 750 1190 1510]; % Selection of
individual frames

f_para.xfact=1000; % Conversion of time unit

f_para.timeunit_string='[s]'; % Time unit, used in labels

f_para.comment_string=['Fig5a']; % Additional comment on the example, will be
used in file and figure names

**Fig. 5B: Further examples of instantaneous
clustering for artificially generated spike trains**

(example 61)

d_para.interval_separators=500:500:d_para.tmax-500;
% Edges
of subsections

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 3 0 0]; % Vector with order of subplots

f_para.plot_mode=2; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.mov_frames=[250 750 1250 1750]; % Selection of individual frames
f_para.xfact=1000; % Conversion of time unit

f_para.timeunit_string='[s]'; % Time unit, used in labels

f_para.comment_string=['Fig5b']; % Additional comment on the example, will be
used in file and figure names

**Fig. 6: Selective temporal averaging**

(example 71)

d_para.interval_separators=500:500:d_para.tmax-500;
% Edges
of subsections

d_para.select_averages={[0 500];[1000 2000];[2000 2500 3000 3500];[0
4000]}; % Selected average over different intervals

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 3 0 0]; % Vector with order of subplots

f_para.plot_mode=2; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.xfact=1000; % Conversion of time unit

f_para.timeunit_string='[s]'; % Time unit, used in labels

f_para.comment_string=['Fig6']; % Additional comment on the example, will be
used in file and figure names

**Fig. 7a: Triggered temporal averaging**

(example 81)

d_para.select_averages={[d_para.tmin d_para.tmax]};
% One continuous intervals for selective temporal
averaging

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 0 0 0]; % Vector with order of subplots

f_para.plot_mode=8; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.dendrograms=1; % Cluster trees from distance
matrices (0-no,1-yes)

f_para.comment_string=['Fig7a']; % Additional comment
on the example, will be used in file and figure names

**Fig. 7b: SPIKE-distance applied to a bivariate
example with varying spike match**

(example 81)

d_para.trigger_averages{1}=spikes(1,1:num_spikes(1));
% Triggered averaging over all time instants when a
certain neuron fires

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 0 0 0]; % Vector with order of subplots

f_para.plot_mode=8; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.dendrograms=1; % Cluster trees from distance
matrices (0-no,1-yes)

f_para.comment_string=['Fig7b']; % Additional comment
on the example, will be used in file and figure names

**Fig. 8: Screenshot from the supplementary movie**

(example 91)

d_para.markers=[500:500:d_para.tmax-500];
% Relevant time
instants

d_para.all_train_group_names={'G1';'G2';'G3';'G4'};
% Names of spike train
groups

d_para.all_train_group_sizes=[10 10 10 10]; % Sizes of respective spike
train groups

d_para.select_averages={500+[0 500 1000 1500]}; % Two
continuous intervals for selective temporal averaging

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 3 0 0]; % Vector with order of subplots

f_para.plot_mode=8; % +1:vs time,+2-different measures and cuts,+4-different
measures,+8-different cuts,+16:different cuts-Movie

f_para.block_matrices=1; % Allows tracing the overall synchronization among
groups of spike trains (0-no,1-yes)

f_para.dendrograms=1; % Cluster trees from distance matrices (0-no,1-yes)

f_para.dendro_color_mode=1; % Coloring of dendrograms
(0-no,1-groups,2-trains)

f_para.comment_string=['Fig8']; % Additional comment on the example, will be
used in file and figure names

**Additional Figure 1: Comparison of piecewise
constant and sampled dissimilarity profile**

(example 21)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2 0 0 0 0]; % Vector with order of subplots

f_para.subplot_size=[1 1.5]; % Vector with relative
size of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.publication=0; % Omits otherwise helpful
information, such as additional comments (0-no,1-yes)

[The dissimilarity
profile with the normal line width is **sampled** (with one value each sample
point), whereas the bold dissimilarity profile is **piecewise constant** (with
one value for each interval in the pooled spike train). This profile is
obtained by representing each interval by the value of its center and
weighting it by its length. This is not only faster, but it actually gives
the exact result, whereas the sampled calculation is a very good
approximation only for sufficiently small sampling intervals (imagine the
example of a rectangular function, at some point any sampled representation
has to cut the right angle). On the other hand, the sampled profile has a
higher temporal resolution in the visualization. If this is not essential or
if the averaged values are sufficient, the calculation of the piecewise
constant distances is preferable.

Reminder:

9-10: ISI-distance (per definition piecewise constant)

11-12: SPIKE-distance (piecewise constant)

13-14: Real-time SPIKE-distance (piecewise constant)

15-16:Future SPIKE-distance (piecewise constant)

17-18: SPIKE-distance (sampled with variable 'time')

19-20: Real-time SPIKE-distance (sampled with variable 'time')

21-22: Future SPIKE-distance (sampled with variable 'time')

In the source codes the sampled variables can be
recognized by the ending '**_t**'.

Regarding the
**plotting parameters**:
The sampled profile was generated because f_para.subplot_posi(18)>0
and the piecewise constant profile was generated because
f_para.subplot_posi(12)>0. Since both values were 2 they were
plotted in the same subplot. And since f_para.subplot_size was set to
[1 1.5], the second subplot is 1.5 times larger than the first.]

**Additional Figure 2: Poisson with different firing
rate - Distance value only depends on distance from diagonal (distances are
scale-free)**

(example 121)

d_para.select_averages={[d_para.tmin d_para.tmax]}; %
Continuous intervals for selective temporal averaging

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 2 0 3 0 0 0 0 0 0 0 0]; % Vector with order of subplots

f_para.plot_mode=8; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.color_norm_mode=3; % normalization of pairwise
color matrices (1-Absolute maximum value,2-Overall maximum value,3-Each one
individually)

**Additional Figure 3: Poisson with different firing
rate - Distance value only depends on distance from diagonal (distances are
scale-free)**

(example 121)

d_para.select_averages={[d_para.tmin d_para.tmax]}; %
Continuous intervals for selective temporal averaging

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 2 0 3 0 0 0 0 0 0 0 0]; % Vector with order of subplots

f_para.plot_mode=8; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.color_norm_mode=3; % normalization of pairwise
color matrices (1-Absolute maximum value,2-Overall maximum value,3-Each one
individually)

**Additional Figure 4: Expectation values for Poisson
spike trains with identical rate (but these values do not depend on the
rate)**

(example 131)

f_para.subplot_posi=[0 1 0 0 0 0 0 0 0 0
0 2 0 3 0 0 0 0 0 0 0 0]; % Vector with order of subplots

f_para.plot_mode=1; % +1:vs time,+2-different measures
and cuts,+4-different measures,+8-different cuts,+16:different cuts-Movie

f_para.publication=0; % Omits otherwise helpful
information, such as additional comments (0-no,1-yes)

[The firing rate of the Poisson spike trains is 1000 Hz so these are almost 1.000.000 spikes. It takes quite some time (~1h on my notebook) to run this example. The interval is cut into smaller segments, and the averaging over all pairs of spike trains is performed for each segment separately.]

See also:

**Measuring spike train synchrony
I: SPIKY (graphical user interface) **

Graphical user interface (Matlab) which can be used to calculate and visualize both the SPIKE- and the ISI-distance between two (or more) spike trains. This is an extension and update of the code in II (there you also find the links to the relevant articles).

**Measuring spike train synchrony
II: cSPIKE **

This Matlab command line software utilizes MEX files with C++ back ends for improved performance and abstraction. It offers the same basic functionality as SPIKY and PySpike and can be easily incorporated as part of your own Matlab program. You can also use cSPIKE as a stand-alone program to calculate the ISI-distance, SPIKE-distance and SPIKE synchronization measures and many more for your spike train sets.

**Measuring spike train synchrony
IV: Event Synchronization**

Matlab-Code to measure the event synchronization and the event delay between two given spike trains

For a detailed description of the method please refer to:

Quian Quiroga R, Kreuz T, and Grassberger P:

Event Synchronization: A simple and fast method to measure synchronicity and time delay patterns.

*Phys.Rev. E* **66**, 041904 (2002) [PDF].

**Measuring spike train synchrony V: Directionality**

Matlab code to calculate the directionality measure L between two given spike trains (or between two continuous datasets or between a spike train and a continuous dataset)

For a detailed description of the method please refer to:

Andrzejak RG, Kreuz T:

Characterizing unidirectional couplings between point processes and flows.

*European Physics Letters* **96**, 50012 (2011) [PDF].

**Measuring spike train synchrony
VI: van Rossum distance and multi-neuron
extension**

Matlab codes to calculate the spike train metric by van Rossum and the multi-neuron extension by Houghton and Sen.

For a detailed description of the method please refer to:

Houghton C, Kreuz T:

On the efficient calculation of van Rossum distances.

*Network: Computation in Neural Systems* **23**, 48 (2012)
[PDF].

**
Measuring spike train synchrony
VII: Victor-Purpura distance
and multi-neuron extension**

Matlab codes to calculate the spike train metric by Victor-Purpura and the multi-neuron extension by Victor-Purpura-Aronov.

(Homepage of Prof. Jonathan D. Victor, Cornell, NY, USA)

See also:

Matlab code for the Victor-Purpura distance which in addition calculates the percentage of spikes that have been matched by a time shift as well as the average time shift

For a detailed description of the algorithm please refer to:

Chicharro D, Kreuz T, Andrzejak RG:

What can spike train distances tell us about the neural code?

*J Neurosci Methods* **199**, 146 (2011)
[PDF].

For **questions and comments** please contact me at "**thomas.kreuz (at) cnr.it**".

These are some **keywords** to help people to find this page.

spike train analysis, measuring spike train synchrony, monitoring spike
train synchrony, ISI-distance, SPIKE-distance, open Matlab source codes,
free download

## Notices

#### FOR SCIENTIFIC USE ONLY

These codes are free of charge for research and education purposes only. Any commercial or military use of this software is prohibited.

The software on this site is provided "as-is," without any expressed or implied warranty. In no event am I or my host institution liable for any damages arising from the use of the software. Since it is distributed for free, I do also not take responsibility for any eventual error in it.