I would like to find bursts among neuronal spikes.

I have an array with neuronal interspike intervals (ISI). I would like to recognise bursts according to the following parameters: a burst starts if the ISI is lower than 6 ms, and ends if the ISI is above 9 ms. The recquirements of a burst is min 3 spikes, and 5 ms in length. There should be 20 ms between 2 bursts.

Can you help me which command should be used to find the bursting spikes? I need the number of bursts.

Answer 1

[MarKf]

There are several features that you may want to take into account depending on your data and on what you want to do. The simplest thing would be to set up a threshold (say, 5 std above, ISI as a diff of spike indexes) that would not capture different types of spikes, tails, and other things, but still that might be enough for you. Have a look at some tools other people have developed for this, for example in fieldtrip (toolbox), here (this has minNspykes, duration and tails) and here (combination of multiple amplitude and spectral features, common approaches).

[MarKf]

Oh sorry nevermind, you said you already have compiled ISI data.

I could whip out something quick (similar to the tools above).

Let's assume sampling freq at 1000 (1ms). Then something like this:

ISIeg = [3 11 12 13 14 40 67 68 ... %example 2 bursts too short (time and spikes)
    99 101 105 114 124 129 134 140 ... % 1 burst, 1 burst too early
    198 224 267 269 276 279 290 301 334 338 345 348 355 ...
    393 450 489 494 530 563]; % then 2 bursts 
bursts = zeros([600, 1]); 
bursts(ISIeg) = 0.5; % to plot bursts & interbursts
minI = 6; maxI = 9; minS = 3; lenB = 5; ibi = 20; %parameters

bibi = ibi*2; % to pad
ISId = diff([-bibi ISIeg ISIeg(end)+bibi]);
startb = ISId < minI;
contb = ISId < maxI;
starts = strfind([false, startb], [0 1]);

Nbursts = 0;
bbepv = -bibi; %initialise previous burst end
for sti = 1:numel(starts)
    bst = starts(sti)-1;
    bis = true; cou = 0;
    while bis
    cou = cou + 1; %next spike
    bend = bst+cou; % cou(nter)+1 how many spikes
    if ~contb(bend+1), bis = false; end % if ISIdiff > maxISI
    end
    bdur = ISIeg(bend)- ISIeg(bst); % burst duration
    bint = ISIeg(bst)-bbepv; % since last burst
    if cou+1 >= minS && bdur >= lenB && bint >= ibi %check "recquirements"
    bbepv = ISIeg(bend);
    bursts(ISIeg(bst):bbepv) = 1; 
    Nbursts = Nbursts+1;
    end
end


plot(bursts), title(['N bursts = ' num2str(Nbursts)]),  
ylim([0 1.5]), ylabel('interbursts = 0.5; bursts = 1')

[Bence Laczó]

Thanks for your anwer, this code looks really awsome. However I don't really understand that in your example ISIeg contains spike positions or interstimulus intervals?

I wanted to attach my ISI data, but it seems that I failed to do it. Now I attached it.

[MarKf]

Actually I see now that the name ISIeg was misleading, sorry, it's not inter-spikes intervals. ISIeg was an example of indices of spikes (when they occur, if you had logical indices instead, i.e. zeros and a 1 every time there is a spike, you could just find them with find(logical_idx) but nevermind that now).

Your shared file is actual ISI probably, so what would become (unpadded) ISId in the previous script. Then it's the same:

ISIdat = load(websave('rd', "https://nl.mathworks.com/matlabcentral/answers/uploaded_files/1317005/ISI.mat"));
ISI = ISIdat.ISI; % this is just load('ISI.mat') file you shared above

spikes = cumsum(ISI); %what was ISIeg
bursts = zeros([sum(ISI), 1]); %or spikes(end)
bursts(spikes) = 0.5; % to plot bursts & interbursts
minI = 6; maxI = 9; minS = 3; lenB = 5; ibi = 20; %parameters

bibi = ibi*2; % to pad
ISId = [bibi ISI' bibi];
startb = ISId < minI;
contb = ISId < maxI;
starts = strfind([false, startb], [0 1]);

Nbursts = 0;
bbepv = -bibi; %initialise previous burst end
for sti = 1:numel(starts)
    bst = starts(sti)-1;
    bis = true; cou = 0;
    while bis
    cou = cou + 1; %next spike
    bend = bst+cou; % cou(nter)+1 how many spikes
    if ~contb(bend+1), bis = false; end % if ISIdiff > maxISI
    end
    bdur = spikes(bend)- spikes(bst); % burst duration
    bint = spikes(bst)-bbepv; % since last burst
    if cou+1 >= minS && bdur >= lenB && bint >= ibi %check "recquirements"
    bbepv = spikes(bend); %in theory you don't have to recreate spikes indices
    bursts(spikes(bst):bbepv) = 1; % you could do everything with ISI
    Nbursts = Nbursts+1; %but the code was already here done, so I guess the above is left as an excercise?
    end
end

plot(bursts), title(['N bursts = ' num2str(Nbursts)]),  
ylim([0 1.5]), ylabel('interbursts = 0.5; bursts = 1')

Make sure the code does what it is supposed to do, with my quick 3 spikes example I created I knew where and how many spikes there were. Here there are many frequent interbusts, so maybe a smaller ibi would make more sense. Consider accepting the answer if that's what you were looking for or if it helped.

[Bence Laczó]

I tried your code but it found 10 bursts among the spikes. With manual counting I found only 4 bursts. I changed the recquirements in the 24th row from if cou+1 >= minS && bdur >= lenB && bint >= ibi to if cou >= minS && bdur >= lenB && bint >= ibi. Now it works fine. Does it make sense? Why did you write cou+1? I don't really understand it.