SpikeTrainSet

Source Code - cSPIKE

Matlab command line library

Description of the codes

cSPIKE v 1.0 29.9.2016: Eero Satuvuori

This is the main part of cSPIKE. The object offers the basic functionality of graphical user interface SPIKY as a command line version. The main functions are implemented with MEX files and C++ backends. Methods are based on the ISI-distance, the SPIKE-distance and the SPIKE-synchronization by Thomas Kreuz et al. The class contains data and methods for storing spike trains as a single set as well as calculating the different distance measures.

To use SpikeTrainSet you need a C++ compiler for your Matlab. Some help

can be found in the following links

https://it.mathworks.com/help/matlab/matlab_external/what-you-need-to-build-mex-files.html

https://it.mathworks.com/support/compilers/R2016a/

SpikeTrainSet objects can be created without compiling MEX files, but

you can't call any distance measures without the appropriate MEX

compiled.

*** Important ***

Before you use cSPIKE for the very first time on a computer you

need to compile the MEX and C++ backends. To do so go to folder

cSPIKEmex and run script MEX_compile.m. This may take few minutes.

Every time you call for MEX files(most functions of the SpikeTrainSet

use them) the folder in which they are located need to be in your

Matlab workspace. Initialize cSPIKE-script will add it to your path. It

needs to be called once every time after starting Matlab and any

additional calls are not needed. This is not called by the SpikeTrainSet!

*****************

The SpikeTrainSet implements the following functions

To create a spike train set object call the constructor:

spiketrains: A cell array with SpikeTrains{1} containing an

array of spike times [spike1 spike2 ...spikeN] for the

first spike train and respectively for the other spike

trains. The object accepts only spike data aligned as row

vectors

beginning: The start time of the recording

endgin: The end time of the recording

STS = SpikeTrainSet(spiketrains, beginning, ending )

Then you can call methods for the spike train set:

time1: Start time of the analysis interval.

time2: End time of the analysis interval. If not given or if

time1 = time2 instantaneous dissimilarity value is given instead.

threshold: Threshold for the adaptive method. If not given the

threshold extracted from the data is used instead.

Profile: Profile object

STS.ISIdistance(time1, time2)

STS.ISIdistanceMatrix(time1, time2)

Profile = STS.ISIdistanceProfile(time1, time2)

STS.AdaptiveISIdistance(time1, time2, threshold)

STS.AdaptiveISIdistanceMatrix(time1, time2,threshold)

Profile = STS.AdaptiveISIdistanceProfile(time1, time2,threshold)

STS.SPIKEdistance(time1, time2)

STS.SPIKEdistanceMatrix(obj, time1, time2)

Profile = STS.SPIKEdistanceProfile(obj, time1, time2)

STS.AdaptiveSPIKEdistance(time1, time2, threshold)

STS.AdaptiveSPIKEdistanceMatrix(obj, time1, time2,threshold)

Profile = STS.AdaptiveSPIKEdistanceProfile(obj, time1, time2,threshold)

STS.RateIndependentSPIKEdistance(time1, time2)

STS.AdaptiveRateIndependentSPIKEdistance(time1, time2, threshold)

STS.SPIKEsynchro(time1,time2)

STS.AdaptiveSPIKEsynchro(time1,time2,threshold)

SPIKESM: Spike synchronization matrix

SPIKEOM: Spike order matrix

normSPIKEOM: normalized spike order matrix

[SPIKESM,SPIKEOM,normSPIKEOM]= STS.SPIKESynchroMatrix(time1, time2)

[SPIKESM,SPIKEOM,normSPIKEOM] = STS.AdaptiveSPIKESynchroMatrix(time1, time2, threshold)

synchro: a cell array identical to spike train set but instead of spike

times it contains SPIKE-synchronization values of each spike

Sorder: a cell array identical to spike train set but instead of spike

times it contains SPIKE-order values of each spike

STOrder: a cell array identical to spike train set but instead of spike

times it contains spike-train-order values of each spike

[synchro,Sorder,STOrder] = STS.SPIKEsynchroProfile(time1, time2)

[synchro,Sorder,STOrder] = STS.AdaptiveSPIKEsynchroProfile(time1, time2,threshold)

varargin: Each time moment is given as a separate parameter. Returns

the dissimilarity averaged over all time points.

STS.TriggeredISImatrix(varargin)

STS.TriggeredAdaptiveISImatrix(threshold, varargin)

STS.TriggeredSPIKEmatrix(varargin)

STS.TriggeredAdaptiveSPIKEmatrix(threshold, varargin)

varargin: Each pair of inputs defines a new interval over which the

distance is defined. A unique union of the intervals is

used when defining the intervals over which the average is

taken.

STS.AveragedISIdistanceMatrix(varargin)

STS.AveragedAdaptiveISIdistanceMatrix(threshold, varargin)

STS.AveragedSPIKEdistanceMatrix(varargin)

STS.AveragedAdaptiveSPIKEdistanceMatrix(threshold, varargin)

Auxiliary functions:

STS.SetData(spiketrains, beginning, ending )

Replaces the data in the object with new data

STS.giveDATA()

Gives the a copy of the data array of the object. Data{1} contains

the edge corrected data and Data{2} the original spikes.

STS.giveTHR()

Returns the threshold value obtained from the set

[time1,time2] = STS.giveTIMES()

Gives the beginning and the end of the recording

STS.plotSpikeTrainSet(colour,widith)

Plots the spike trains to current axis with colour and spike widith

given. Default colour is black.

BSD license:

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the author nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL HOLDER BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Class for enclosing spiketrains into one class and handling the set.
PRIVATE METHODS **************************************************
Checking data availability
Unique Union
Checking that the DATA given is valid
Cutting function
Check input times
Correct the DATA Edges
Form pooled spiketrain of all spiketrains
Time shift functions.
PUBLIC METHODS ***************************************************
Constructor
Setting data
ISI METHODS ----------------------------------------------------
ISI-distance for interval/point
Adaptive ISI-distance for interval/point
ISI-distance/dissimilarity matrix
Adaptive ISI-distance/dissimilarity matrix
Triggered ISI-dissimilarity matrix
Triggered Adaptive ISI-dissimilarity matrix
Averaged ISI-distance matrix
Averaged Adaptive ISI-distance matrix
ISI-distance profile
Adaptive ISI-distance profile
SPIKE METHODS --------------------------------------------------
SPIKE-distance for interval/point
Adaptive SPIKE-distance for interval/point
Rate independent SPIKE distance for interval/point
Adaptive Rate independent SPIKE distance for interval/point
SPIKE distance/dissimilarity matrix
Adaptive SPIKE distance/dissimilarity matrix
Triggered SPIKE-dissimilarity matrix
Triggered Adaptive SPIKE-dissimilarity matrix
Averaged SPIKE-distance matrix
Averaged Adaptive SPIKE-distance matrix
SPIKE-distance profile
Adaptive SPIKE-distance profile
SPIKE-Synchronization METHODS ----------------------------------
SPIKE-Synchronization
Adaptive SPIKE-synchronization
SPIKE-synchronization profile
Adaptive SPIKE-synchronization profile
SPIKE-synchronization matrix
Adaptive SPIKE-synchronization matrix
Auxiliary functions --------------------------------------------
Get function
Threshold value extracted from the data
Start and end of the data set
Plotting function

Class for enclosing spiketrains into one class and handling the set.

classdef SpikeTrainSet < handle

properties (Access = private)

% Internal datastores.

DataArray;

PooledArray;

timeStart;

timeEnd;

threshold;

timeShift;

realISIs;

% Internal accuracy used by the program. This is the maximum

% accuracy used. Increasing the accuracy may lead to malfunctions.

accuracy = 10^-12;

end

PRIVATE METHODS **************************************************

methods (Access = private)

Checking data availability

Function is given a bool value if the information is crucial for the calling function. In case it is an error is thrown. If it is not, just a bool value of 0 is returned if there is no data yet.

function bool = HaveData(obj, crucial)

% Returning 1 if there is data and 0 if there is not.

if isempty(obj.DataArray)

bool = 0;

if crucial

error('No data available');

end

else

bool = 1;

end

Unique Union

Takes a set of intervals as a cell array so that each pair of two values is an interval. The output is a unique union of the intervals.

function ReturnArray = UniqueUnion(obj, CellArray )

starts = [];

ends = [];

for i = 1:2:length(CellArray)

starts(end+1) = CellArray{i};

ends(end+1) = CellArray{i+1};

end

UniqueArray = 0;

while ~UniqueArray

moves = 0;

% Taking each first spike and if it is in an interval, move

% it to the border of the interval.

for i = 1:length(starts)

for ii = 1:length(starts)

if (i ~= ii)

if (starts(i) > starts(ii) && starts(i)< ends(ii))

starts(i) = starts(ii);

moves = moves + 1;

end

% Taking each last spike and doing the same.

for i = 1:length(ends)

for ii = 1:length(ends)

if (i ~= ii)

if (ends(i) > starts(ii) && ends(i) < ends(ii))

ends(i) = ends(ii);

moves = moves + 1;

end

% If there were no moves this round, all the spikes are at

% the edges of their intervals. There may be multiple same

% intervals.

if moves == 0

% now the array is unique. No need to continue.

UniqueArray = 1;

% Taking each interval only once.

Ustarts = unique(starts);

Uends = unique(ends);

% Since the intervals are in time order and unique,

% they will be in correct time order after unique

% command. The output is formed based on the intervals

ReturnArray = cell(0);

for i = 1:length(Ustarts)

ReturnArray{end+1} = Ustarts(i);

ReturnArray{end+1} = Uends(i);

end

Checking that the DATA given is valid

Function checks the data and gives an error if the data is not in a correct form

function CheckData(obj,Data)

if ~iscell(Data)

error('Data is not a cell array');

end

if length(Data) == 1

error('One spiketrain');

end

for i = 1:length(Data)

if ndims(Data)~=2 || size(Data,2) ~= 1;

error('Data has incorrectly aligned cell array');

end

if ndims(Data{i}) ~=2 && (size(Data{i},1)||size(Data{i},2))

error('Data has too many dimensions');

end

if size(Data{i},1) ~= 1

if ~isempty(Data{i})

error('Data has incorrectly aligned vectors');

end

if ~isnumeric(Data{i})

error('Data is not all numeric');

end

Cutting function

This function cuts the data so that spikes outside of boundaries are removed.

function Data = CutData(obj,Data, START, END)

for i = 1:length(Data)

inds = not(abs(sign(sign(START-obj.accuracy - Data{i})...

+ sign(END+obj.accuracy - Data{i}))));

Data{i} = Data{i}(inds);

end

function Data = RoundData(obj, Data)

dt = obj.accuracy;

for i = 1:length(Data)

Data{i} = round(Data{i}/dt)*dt;

end

Check input times

Checking that times given are in correct order and within boundaries

function CheckTime(obj,Time1,Time2)

if Time1 > Time2

error('Time interval boundaries incorrect');

end

if ~((Time1 >= obj.timeStart) && (Time2 <= obj.timeEnd))

error('time interval not within data bounds');

end

Correct the DATA Edges

Adds edge correction and makes sure that there are only 1 spike at each time instant. Also rearranges the data into time series.

function Corrected = EdgeCorrection(obj,Data,time1,time2)

PooledISIs = [];

Corrected = cell(1,length(Data))';

for i = 1:length(Data)

% 0 or 1 spikes cases

if isempty(Data{i}) || length(Data{i}) == 1

if isempty(Data{i})

first = -(time2-time1)*1.2;

last = (time2-time1)*2.2;

Corrected{i} = [first obj.timeStart obj.timeEnd last];

Corrected{i} = unique(Corrected{i});

thisST = obj.timeEnd-obj.timeStart;

else

fullST = 0; %flag

first = obj.timeStart;

last = obj.timeEnd;

if Data{i} == first

first = -(time2-time1)*1.2;

thisST = obj.timeEnd-obj.timeStart;

fullST = 1;

end

if Data{i} == last

last = (time2-time1)*2.2;

thisST = obj.timeEnd-obj.timeStart;

fullST = 1;

end

if ~fullST

thisST = [last-Data{i} Data{i}-first];

end

Corrected{i} = [first Data{i} last];

Corrected{i} = unique(Corrected{i});

end

else

%The rest

thisST = Data{i}(2:end) - Data{i}(1:end-1);

if(Data{i}(1) ~= time1)

first = Data{i}(1) - max( Data{i}(2)-Data{i}(1), Data{i}(1) - obj.timeStart );

thisST = [thisST Data{i}(1)-first];

else

% In case the first spike is at the border it is

% not treated as auxiliary spike. For calculation

% the next spike does not matter so we add an

% "auxiliary" spike very far from the dataset to

% mark that the border is a real one.

first = -(time2-time1)*1.2;

end

if(Data{i}(end)~= time2)

last = Data{i}(end) + max( Data{i}(end)-Data{i}(end-1),obj.timeEnd - Data{i}(end));

thisST = [thisST last-Data{i}(end)];

else

% Same for the other end.

last = (time2-time1)*2.2;

end

Corrected{i} = [first last Data{i}];

% At the end making sure that each dataset has only one

% spike at each time instant and that they are in time

% order

Corrected{i} = unique(Corrected{i});

end

PooledISIs = [PooledISIs thisST];

end

if(isempty (PooledISIs))

obj.threshold = 0;

else

obj.threshold = sqrt(mean(PooledISIs.^2));

end

obj.realISIs = PooledISIs;

end

Form pooled spiketrain of all spiketrains

Pools all spiketrains into one. Uses unique, so that multiple spikes at same times are filtered out.

function Pooled = Pool(obj , Data)

Pooled = [];

for i = 1:length(Data)

Pooled = [Pooled Data{i}];

end

Pooled = unique(Pooled);

end

Time shift functions.

These functions take care of the difference between "real" times and internal times. The object computes using positive spiketimes and in case the recording begins at negative times, all data is shifted so that the recording begins from 0 internally. For output the times are again shifted back to original.

function [time1,time2,Data] = CreateTimeShift(obj,time1,time2,Data)

obj.timeShift = time1;

time1 = time1-obj.timeShift;

time2 = time2-obj.timeShift;

for i = 1:length(Data)

Data{i}=Data{i}-obj.timeShift;

end

function [time] = InputTimeShift(obj,time)

time = time-obj.timeShift;

end

function [time] = outputTimeShift(obj,time)

time = time+obj.timeShift;

end

PUBLIC METHODS ***************************************************

methods (Access = public)

Constructor

Builds a spiketrain set with the spiketrains given and beginning and end times.

function obj = SpikeTrainSet(spiketrains, beginning, ending )

obj.timeShift = 0;

obj.SetData(spiketrains, beginning, ending );

end

Not enough input arguments.

Error in SpikeTrainSet (line 471)

obj.SetData(spiketrains, beginning, ending );

Setting data

Function sets data to the classes own variables. Calls appropriate checking functions to input. Data is only accepted in a single cell array containing spiketrains. If called to a spiketrain set that has already been set, the user is asked for input if it is to be replaced.

function SetData(obj, Data, time1, time2)

Data = Data';

% checking that the data is in correct form.

obj.CheckData(Data);

%Rounding data to 12 digits accuracy

Data = obj.RoundData(Data);

% Cutting spikes outside of interval

Data = obj.CutData(Data, time1, time2);

% If there was no data, add the data

if obj.HaveData(0)

%presetting answer variable into not n or y

answer = 'h';

% If there was data, ask if the data is to be replaced

while answer ~= 'y' && answer ~= 'Y' && answer ~= 'n' && answer ~= 'N'

answer = input('Do you want to replace the data?(y/n)','s');

end

%

if answer == 'y' || answer ~= 'Y'

% Swiping old arrays clean

obj.DataArray = [];

obj.PooledArray= [];

% Adding new values to the array

obj.DataArray{2} = Data;

[time1,time2,Data] = obj.CreateTimeShift(time1,time2,Data);

obj.timeStart = time1;

obj.timeEnd = time2;

obj.DataArray{1} = obj.EdgeCorrection(Data);

obj.PooledArray = obj.Pool(Data);

end

else

% Setting the arrays for the first time

obj.DataArray{2} = Data;

[time1,time2,Data] = obj.CreateTimeShift(time1,time2,Data);

obj.timeStart = time1;

obj.timeEnd = time2;

obj.DataArray{1} = obj.EdgeCorrection(Data,time1,time2);

obj.PooledArray = obj.Pool(Data);

end

ISI METHODS ----------------------------------------------------

The methods in this section are applying ISI distance measure to the data. The section also contains Adaptive extension.

ISI-distance for interval/point

Calculates ISI distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, return dissimilarity of the spiketrains at that time instant.

function ISId = ISIdistance(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

ISId = AdaptiveISIdistance(obj, time1, time2, 0);

end

Adaptive ISI-distance for interval/point

Calculates SPIKE distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, return dissimilarity of the spiketrains at that time instant.

function ISId = AdaptiveISIdistance(obj, time1, time2, threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin == 2

time2 = time1;

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

if(time1 == time2)

ISId = mexAdaptiveISIDistance(obj.DataArray{1},threshold,time1);

else

ISId = mexAdaptiveISIDistance(obj.DataArray{1},threshold,time1, time2);

end

ISI-distance/dissimilarity matrix

Returns the distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs.

function ISIdM = ISIdistanceMatrix(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

ISIdM = AdaptiveISIdistanceMatrix(obj, time1, time2, 0);

end

Adaptive ISI-distance/dissimilarity matrix

Returns the adaptive distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs. Threshold can be given as a parameter. If not given the one extracted from the data is used.

function ISIdM = AdaptiveISIdistanceMatrix(obj, time1, time2,threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin == 2

time2 = time1;

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

ISIdM = mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold ,time1, time2);

end

Triggered ISI-dissimilarity matrix

Returns the averaged dissimilarity matrix of the spike trains averaged over time instants defined in inputs. Each instant is given as new variables: (1,2,5,6) being average over time instants 1,2,5 and 6.

function ISIdM = TriggeredISImatrix(obj, varargin)

if nargin < 3

error('Not enough input arguments');

end

ISIdM = TriggeredAdaptiveISImatrix(obj,0, varargin);

end

Triggered Adaptive ISI-dissimilarity matrix

Returns the averaged dissimilarity matrix of the spike trains averaged over time instants defined in inputs. Each instant is given as new variables: (1,2,5,6) being average over time instants 1,2,5 and 6. If a negative threshold value is given the data exctracted threshold is used.

function ISIdM = TriggeredAdaptiveISImatrix(obj,threshold, varargin)

if threshold < 0

threshold = obj.threshold;

end

if nargin < 4

error('Not enough input arguments');

end

obj.HaveData(1);

time = obj.InputTimeShift(varargin{1});

obj.CheckTime(time,time);

ISIdM = mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold,time, time);

for i = 2:(nargin-1)

time = obj.InputTimeShift(varargin{i});

obj.CheckTime(time,time);

ISIdM = ISIdM + mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold,time, time);

end

ISIdM = ISIdM./(nargin-2);

end

Averaged ISI-distance matrix

Returns the distance matrix of the spike trains averaged over time intervals defined as pairs of inputs. The inputs must be in pairs. If the intervals overlap, the average is computed over unique union of the intervals. Input is given as new pairs of variables: (1,2,5,6) being average over time intervals [1,2] and [5,6].

function ISIdM = AveragedISIdistanceMatrix(obj, varargin)

if nargin < 3

error('Not enough input arguments');

end

if mod(nargin,2) == 0;

error('Not a paired number of time points');

end

ISIdM = AveragedAdaptiveISIdistanceMatrix(obj,0, varargin);

end

Averaged Adaptive ISI-distance matrix

Returns the distance matrix of the spike trains averaged over time intervals defined as pairs of inputs. The inputs must be in pairs. If the intervals overlap, the average is computed over unique union of the intervals. Input is given as new pairs of variables: (1,2,5,6) being average over time intervals [1,2] and [5,6]. Threshold can be given as a parameter. If a negative threshold value is given the data exctracted threshold is used.

function ISIdM = AveragedAdaptiveISIdistanceMatrix(obj,threshold, varargin)

if threshold < 0

threshold = obj.threshold;

end

if nargin < 4

error('Not enough input arguments');

end

if mod(nargin,2) == 1;

error('Not a paired number of time points');

end

obj.HaveData(1);

% Forming unique union of input

Timeslices = obj.UniqueUnion(varargin);

time1 = Timeslices{1};

time2 = Timeslices{2};

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

obj.CheckTime(time1,time2);

TotalTime = time2-time1;

ISIdM = (time2-time1)*mexISIistanceMatrix(obj.DataArray{1},threshold,time1, time2);

for i = 3:2:length(Timeslices)

time1 = Timeslices{i};

time2 = Timeslices{i+1};

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

obj.CheckTime(time1,time2);

ISIdM = ISIdM + (time2-time1)*mexISIistanceMatrix(obj.DataArray{1},threshold,time1, time2);

TotalTime = TotalTime + (time2-time1);

end

% averaging over time. Each matrix is weighted

ISIdM = ISIdM./TotalTime;

end

ISI-distance profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting. See object Profile for more information.

function ISIp = ISIdistanceProfile(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

ISIp = AdaptiveISIdistanceProfile(obj, time1, time2,0);

end

Adaptive ISI-distance profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting. See object Profile for more information.

function ISIp = AdaptiveISIdistanceProfile(obj, time1, time2,threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

internalProfile = mexAdaptiveISIDistanceProfile(obj.DataArray{1},threshold,time1, time2);

internalProfile(2,:) = obj.outputTimeShift(internalProfile(2,:));

ISIp = Profile(internalProfile);

end

SPIKE METHODS --------------------------------------------------

The methods in this section are applying the SPIKE-distance measure to the data. The section also contains Adaptive and rate intedependent extensions.

SPIKE-distance for interval/point

Calculates SPIKE distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, return dissimilarity of the spiketrains at that time instant.

function SPIKEd = SPIKEdistance(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

SPIKEd = AdaptiveSPIKEdistance(obj, time1, time2, 0);

end

Adaptive SPIKE-distance for interval/point

Calculates SPIKE distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, return dissimilarity of the spiketrains at that time instant.

function SPIKEd = AdaptiveSPIKEdistance(obj, time1, time2, threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin == 2

time2 = time1;

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

if(time1 == time2)

SPIKEd = mexAdaptiveSPIKEDistance(obj.DataArray{1},threshold,time1);

else

SPIKEd = mexAdaptiveSPIKEDistance(obj.DataArray{1},threshold,time1, time2);

end

Rate independent SPIKE distance for interval/point

Calculates SPIKE distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, returns dissimilarity of the spiketrains at that time instant.

function SPIKEd = RateIndependentSPIKEdistance(obj, time1, time2)

SPIKEd = AdaptiveRateIndependentSPIKEdistance(obj, time1, time2, 0);

end

Adaptive Rate independent SPIKE distance for interval/point

Calculates SPIKE distance for given interval between time1 and time2. time1 < time2. If no input times are given return value for whole data set interval that is defined. If only one time is given or both are the same, returns dissimilarity of the spiketrains at that time instant. This will use the predetermined threshold. If you then want to set threshold yourself give it as third parameter.

function SPIKEd = AdaptiveRateIndependentSPIKEdistance(obj, time1, time2, threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

if(time1 == time2)

SPIKEd = mexAdaptiveRateIndependentSPIKEDistance(obj.DataArray{1},threshold,time1);

else

SPIKEd = mexAdaptiveRateIndependentSPIKEDistance(obj.DataArray{1},threshold,time1, time2);

end

SPIKE distance/dissimilarity matrix

Returns the distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs.

function SPIKEdM = SPIKEdistanceMatrix(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

SPIKEdM = AdaptiveSPIKEdistanceMatrix(obj, time1, time2, 0);

end

Adaptive SPIKE distance/dissimilarity matrix

Returns the distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs.

function SPIKEdM = AdaptiveSPIKEdistanceMatrix(obj, time1, time2,threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin == 2

time2 = time1;

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

SPIKEdM = mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);

end

Triggered SPIKE-dissimilarity matrix

Returns the averaged dissimilarity matrix of the spike trains averaged over time instants defined in inputs. Each instant is given as new variables: (1,2,5,6) being average over time instants 1,2,5 and 6.

function SPIKEdM = TriggeredSPIKEmatrix(obj, varargin)

if nargin < 3

error('Not enough input arguments');

end

SPIKEdM = TriggeredAdaptiveSPIKEmatrix(obj,0, varargin);

end

Triggered Adaptive SPIKE-dissimilarity matrix

Returns the averaged dissimilarity matrix of the spike trains averaged over time instants defined in inputs. Each instant is given as new variables: (1,2,5,6) being average over time instants 1,2,5 and 6. If a negative threshold value is given the data exctracted threshold is used.

function SPIKEdM = TriggeredAdaptiveSPIKEmatrix(obj,threshold, varargin)

if threshold < 0

threshold = obj.threshold;

end

if nargin < 4

error('Not enough input arguments');

end

obj.HaveData(1);

time = obj.InputTimeShift(varargin{1});

obj.CheckTime(time,time);

SPIKEdM = mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time, time);

for i = 2:(nargin-1)

time = obj.InputTimeShift(varargin{i});

obj.CheckTime(time,time);

SPIKEdM = SPIKEdM + mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time, time);

end

SPIKEdM = SPIKEdM./(nargin-1);

end

Averaged SPIKE-distance matrix

Returns the distance matrix of the spike trains averaged over time intervals defined as pairs of inputs. The inputs must be in pairs. If the intervals overlap, the average is computed over unique union of the intervals. Input is given as new pairs of variables: (1,2,5,6) being average over time intervals [1,2] and [5,6].

function SPIKEdM = AveragedSPIKEdistanceMatrix(obj, varargin)

if nargin < 3

error('Not enough input arguments');

end

if mod(nargin,2) == 0;

error('Not a paired number of time points');

end

SPIKEdM = AveragedAdaptiveSPIKEdistanceMatrix(obj, 0, varargin);

end

Averaged Adaptive SPIKE-distance matrix

Returns the distance matrix of the spike trains averaged over time intervals defined as pairs of inputs. The inputs must be in pairs. If the intervals overlap, the average is computed over unique union of the intervals. Input is given as new pairs of variables: (1,2,5,6) being average over time intervals [1,2] and [5,6]. Threshold can be given as a parameter. If a negative threshold value is given the data exctracted threshold is used.

function SPIKEdM = AveragedAdaptiveSPIKEdistanceMatrix(obj, threshold, varargin)

if threshold < 0

threshold = obj.threshold;

end

if nargin < 4

error('Not enough input arguments');

end

if mod(nargin,2) == 1;

error('Not a paired number of time points');

end

obj.HaveData(1);

% Forming unique union of input

Timeslices = obj.UniqueUnion(varargin);

time1 = Timeslices{1};

time2 = Timeslices{2};

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

obj.CheckTime(time1,time2);

TotalTime = time2-time1;

SPIKEdM = (time2-time1)*mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);

for i = 3:2:length(Timeslices)

time1 = Timeslices{i};

time2 = Timeslices{i+1};

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

obj.CheckTime(time1,time2);

SPIKEdM = SPIKEdM + (time2-time1)*mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);

TotalTime = TotalTime + (time2-time1);

end

% averaging over time. Each matrix is weighted

SPIKEdM = SPIKEdM./TotalTime;

end

SPIKE-distance profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting.

function SPIKEp = SPIKEdistanceProfile(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

SPIKEp = AdaptiveSPIKEdistanceProfile(obj, time1, time2, 0);

end

Adaptive SPIKE-distance profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting.

function SPIKEp = AdaptiveSPIKEdistanceProfile(obj, time1, time2,threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

internalProfile = mexAdaptiveSPIKEDistanceProfile(obj.DataArray{1},threshold,time1, time2);

internalProfile(2,:) = obj.outputTimeShift(internalProfile(2,:));

SPIKEp = Profile(internalProfile);

end

SPIKE-Synchronization METHODS ----------------------------------

The methods in this section are applying SPIKE-synchronization measure to the data

SPIKE-Synchronization

This function returns the value of the method from the interval given in the input.

function synchro = SPIKEsynchro(obj,time1,time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

synchro = AdaptiveSPIKEsynchro(obj,time1,time2, 0);

end

Adaptive SPIKE-synchronization

This function returns the value of the method from the interval given in the input. If no threshold value is given the data exctracted value is used.

function synchro = AdaptiveSPIKEsynchro(obj,time1,time2,threshold)

obj.HaveData(1);

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

synchro = mexAdaptiveSPIKESynchro(obj.DataArray{1},threshold/2,obj.timeStart,obj.timeEnd,time1,time2);

end

SPIKE-synchronization profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting.

function [synchro, Sorder, STOrder] = SPIKEsynchroProfile(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

[synchro, Sorder, STOrder] = AdaptiveSPIKEsynchroProfile(obj, time1, time2, 0);

end

Adaptive SPIKE-synchronization profile

Forms a profile between times 1 and 2. The output is a Matlab class Profile, which contains the profile and can even perform a simple plotting.

function [synchro, Sorder, STOrder] = AdaptiveSPIKEsynchroProfile(obj, time1, time2,threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

% Change later into a real profile

[synchro, Sorder, STOrder] = mexAdaptiveSPIKESynchroProfile(obj.DataArray{1},threshold/2,obj.timeStart,obj.timeEnd,time1, time2);

end

SPIKE-synchronization matrix

Returns the distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs.

function [SPIKESM,SPIKEOM,normSPIKEOM] = SPIKESynchroMatrix(obj, time1, time2)

if nargin == 1

time1 = outputTimeShift(obj,obj.timeStart);

time2 = outputTimeShift(obj,obj.timeEnd);

end

if nargin == 2

time2 = time1;

end

[SPIKESM,SPIKEOM,normSPIKEOM] = AdaptiveSPIKESynchroMatrix(obj, time1, time2, 0);

end

Adaptive SPIKE-synchronization matrix

Returns the distance matrix of the spike trains averaged over time interval defined by time 1 and time 2. If not given, averaged over whole interval. If given only one time or the same one twice, the returned matrix is dissimilarity between all pairs. If no threshold value is given thedata exctracted value is used.

function [SPIKESM,SPIKEOM,normSPIKEOM] = AdaptiveSPIKESynchroMatrix(obj, time1, time2, threshold)

if nargin == 1

time1 = obj.timeStart;

time2 = obj.timeEnd;

else

time1 = obj.InputTimeShift(time1);

time2 = obj.InputTimeShift(time2);

end

if nargin == 2

time2 = time1;

end

if nargin < 4

threshold = obj.threshold;

end

obj.HaveData(1);

obj.CheckTime(time1,time2);

NroSpiketrains = length(obj.DataArray{1});

SPIKESM = ones(NroSpiketrains);

SPIKEOM = zeros(NroSpiketrains);

normSPIKEOM = zeros(NroSpiketrains);

for i = 1:NroSpiketrains-1

for ii = i:NroSpiketrains

if i ~= ii

pair{1} = obj.DataArray{2}{i};

pair{2} = obj.DataArray{2}{ii};

spikes = length(obj.DataArray{1}{i})+length(obj.DataArray{1}{ii})-4;

if spikes > 0

STS = SpikeTrainSet(pair,obj.timeStart,obj.timeEnd);

SPIKESM(i,ii) = STS.AdaptiveSPIKEsynchro(time1, time2, threshold);

SPIKESM(ii,i) = SPIKESM(i,ii);

[~,~,profile] = STS.AdaptiveSPIKEsynchroProfile(time1, time2, threshold);

SPIKEOM(i,ii) = sum(profile{1});

SPIKEOM(ii,i) = sum(profile{2});

normSPIKEOM(i,ii) = SPIKEOM(i,ii)/spikes;

normSPIKEOM(ii,i) = SPIKEOM(ii,i)/spikes;

end

Auxiliary functions --------------------------------------------

Get function

Gives a copy of the spike data contained in the spike train set

function DATA = giveDATA(obj)

obj.HaveData(1);

DATA = obj.DataArray;

end

Threshold value extracted from the data

returns the threshold value extracted from this dataset

function [THR, Mean] = giveTHR(obj)

obj.HaveData(1);

THR = obj.threshold;

Mean = mean(obj.realISIs);

end

Start and end of the data set

returns an array with start and end times of the recording.

function [time1,time2] = giveTIMES(obj)

obj.HaveData(1);

time1 = obj.outputTimeShift(obj.timeStart);

time2 = obj.outputTimeShift(obj.timeEnd);

end

Plotting function

Plots the spike trains to current axis

function plotSpikeTrainSet(obj,colour,width)

obj.HaveData(1);

if nargin == 1

colour = 'k';

end

if nargin < 3

width = 1;

end

beginning = obj.outputTimeShift(obj.timeStart);

ending = obj.outputTimeShift(obj.timeEnd);

hold on;

for i = 1:length(obj.DataArray{2})

SpikeTrain = obj.DataArray{2}{i};

plotSpikeTrain(SpikeTrain,0.8,(length(obj.DataArray{2})+1-i)-0.5,colour,width)

end

axis([beginning,ending,0,length(obj.DataArray{2})])

set(gca,'ytick',[]);

end

Published with MATLAB® R2015b

SpikeTrainSet

Source Code - cSPIKE

Matlab command line library

Description of the codes

cSPIKE v 1.0 29.9.2016: Eero Satuvuori

Contents

Class for enclosing spiketrains into one class and handling the set.

PRIVATE METHODS **************************************************

Checking data availability

Unique Union

Checking that the DATA given is valid

Cutting function

Check input times

Correct the DATA Edges

Form pooled spiketrain of all spiketrains

Time shift functions.

PUBLIC METHODS ***************************************************

Constructor

Setting data

ISI METHODS ----------------------------------------------------

ISI-distance for interval/point

Adaptive ISI-distance for interval/point

ISI-distance/dissimilarity matrix

Adaptive ISI-distance/dissimilarity matrix

Triggered ISI-dissimilarity matrix

Triggered Adaptive ISI-dissimilarity matrix

Averaged ISI-distance matrix

Averaged Adaptive ISI-distance matrix

ISI-distance profile

Adaptive ISI-distance profile

SPIKE METHODS --------------------------------------------------

SPIKE-distance for interval/point

Adaptive SPIKE-distance for interval/point

Rate independent SPIKE distance for interval/point

Adaptive Rate independent SPIKE distance for interval/point

SPIKE distance/dissimilarity matrix

Adaptive SPIKE distance/dissimilarity matrix

Triggered SPIKE-dissimilarity matrix

Triggered Adaptive SPIKE-dissimilarity matrix

Averaged SPIKE-distance matrix

Averaged Adaptive SPIKE-distance matrix

SPIKE-distance profile

Adaptive SPIKE-distance profile

SPIKE-Synchronization METHODS ----------------------------------

SPIKE-Synchronization

Adaptive SPIKE-synchronization

SPIKE-synchronization profile

Adaptive SPIKE-synchronization profile

SPIKE-synchronization matrix

Adaptive SPIKE-synchronization matrix

Auxiliary functions --------------------------------------------

Get function

Threshold value extracted from the data

Start and end of the data set

Plotting function