Source Code - Measure of directionality L



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.
Eur Phys Lett 96, 50012 (2011) [PDF].

[Andrzejak_Kreuz_Source_Code.zip]          Copyright:  Ralph G. Andrzejak, Thomas Kreuz

This zip-file includes all of these files:

[AndrzejakKreuzExample.m]

[AndrzejakKreuzL.m]

[AndrzejakKreuzDflow.m]

[AndrzejakKreuzDpointprocess.m]

[f_SPIKE_ISI_distance_new.m]

[AndrzejakKreuzExampleData.mat]

In order to understand the following documentation, you should at first read the paper. In case you have questions, please contact ralphandrzejak(at)yahoo.de (see also his webpage).

Below we use the example of the coupled Hindmarsh-Rose dynamics studied in the paper to illustrate the code. In particular we have:

Signal1: A flow of the dynamics X

Signal2: A flow of a replica surrogate of the dynamics X

Signal3: A point process derived from dynamics Y.

All signals are of length 200T (see again the paper for the definition of T). The dynamics X is driving dynamics Y with epsilon = 0.1442. The notation for the following source code parameters corresponds to the one used in the paper.

q: window length

s: step size

W: Theiler window

Q: total duration,

N: Number of windows

k: Number of nearest neighbours

delta_t: sampling step

Remark on time unit convention: The time step between subsequent entries in the Matlab arrays of flow signals corresponds to delta_t. In other words, by going from the element with index i to the element with index i+1, we make a step of delta_t.

Accordingly, the parameters q, s, and Q are specified in units of delta_t. The parameters q, s and multiples thereof can directly be used as array indices, and Q is the length of the flow arrays. Furthermore, the spike times of the point process signals are also specified in units of delta_t and therefore one should use delta_t=1 as input for the function ‘AndrzejakKreuzDpointprocess’ below. Accordingly, the length of the flow signals is 400.000 samples. And the spike times are bounded by 0 and 400.000. N and W are in units of time windows.

To run the code copy: AndrzejakKreuzExampleData.mat,

AndrzejakKreuzExample.m, AndrzejakKreuzL.m, AndrzejakKreuzDflow.m,

AndrzejakKreuzDpointprocess.m, f_SPIKE_ISI_distance_new.m into some directory and call AndrzejakKreuzExample.m Please look through the code and read the following short documentation of the main functions.

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

 A distance matrix from a flow is calculated using:

function DF=AndrzejakKreuzDflow(flow,Q,q,s,W)

Here flow is a 1-D vector of length Q containing the sample values of the flow. The output DF is a (N,N) distance matrix. All other parameters are explained above.

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

A distance matrix from a point process is calculated using

function DP=AndrzejakKreuzDpointprocess(spiketimes,Q,q,s,W,dt,dchoice)

spiketimes is a 1-D vector of containing the times of spikes. The length of the vector is determined by the number of spikes.

dchoice: In the paper we used the ISI distance. With this code you can also use the

SPIKE distance (see references below for more information on these distances.).

Depending on the value of dchoice, you get as output the ISI-distance, the SPIKE- distance or both:

If you set

dchoice=1: DP is a 2-D array of size (N,N) containing the SPIKE distance.

dchoice=2: DP is a 2-D array of size (N,N) containing the ISI distance.

dchoice=2: DP is a 3-D array of size (2,N,N) containing the SPIKE distance in

DP(1,:,:) and the ISI distance in DP(2,:,:).

All other parameters are explained above.

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

The main function is:

function L = AndrzejakKreuzL(DArray,k,W)

DArray is a 3-D array of size (N,N,ND). It contains ND distance matrices of size (N,N).

In our example we have 3 signals. Accordingly, ND=3. For your own data you can use any number of ND>1 distance matrices (At some point limited by the memory of your computer). You do not have to specify ND. It will be detected automatically from the

size of DArray. The individual distance matrices can be obtained from flows and/or point processes. The order used in our example is arbitrary. Whatever order you choose, it will also be used for the output:

L: a 2-D array of size (ND,ND). In element L(a,b)we have L(a|b). For example L(2,3)contains L(signal2 | signal3) = L(replica surrogate of flow X | point process derived from dynamics Y). The diagonal is set to zero. For the present example, since X is driving Y, we get a high value for L(signal 1 | signal 3) = L(flow X | point process derived from dynamics Y). Lower values are obtained for the other non-diagonal entries.

 

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).

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). Also submitted to the arXiv [PDF].

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].

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 III:   SPIKE- and ISI-Distance

Matlab codes to calculate both the SPIKE- and the 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].

Also submitted to the arXiv [PDF].

Kreuz T:
SPIKE-distance.
Scholarpedia
7(12), 30652 (2012).

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 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".


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 I am 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.