Publish code, license and documention

The files are related to the publication Välkki IA, Lenk K, Mikkonen JE, Kapucu FE and Hyttinen JAK (2017) Network-Wide Adaptive Burst Detection Depicts Neuronal Activity with Improved Accuracy. Front. Comput. Neurosci. 11:40. doi: 10.3389/fncom.2017.00040

CalculateCMANoRound.m

@@ -0,0 +1,80 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% CalculateCMANoRound.m: Runs the CMA burst detection
+% authors: Inkeri A. Välkki, Kerstin Lenk, Jarno E. Mikkonen, Fikret E. Kapucu, Jari A. K. Hyttinen
+% date: 2016 - 2019
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [BurstInfo,CMAinfo, ISIinfo] = CalculateCMANoRound(CMAtype, SpikeTimes, IDs, ST_Duration_ms)
+% [BurstInfo, CMAinfo, ISIinfo] = CalculateCMANoRound( CMAtype, SpikeTrains, IDs, ST_Duration_ms)
+% in:
+%   CMAtype: 'original', 'net_alpha' or 'net'
+%   SpikeTimes: spike times in ms, a cell with individual spike sequences as
+%   vectors.
+%   IDs: some IDs for the spike trains (numbers, channels...)
+%   ST_Duration_ms: Length of the spike sequences in ms, one per spike
+%   sequence or one common for all
+% out:
+%   BurstInfo: bursting statistics as a table
+%   CMAinfo: the parameters that were used in burst detection as table
+%   ISIinfo: statistics of the ISIs
+
+%%
+
+% ST length in min
+ST_Duration_min = ST_Duration_ms / 1000 / 60;
+
+% Get ISI and CMA infos needed
+[CMAinfo, ISIinfo] = getCMAinfo(CMAtype, SpikeTimes, IDs);
+
+% Go through all spike sequences
+for s = 1:length(SpikeTimes)
+
+    % Spike Sequence to analyze
+    Spikes = SpikeTimes{s};
+
+    if length(ST_Duration_min) > 1
+        simtime = ST_Duration_min(s);
+    else
+        simtime = ST_Duration_min;
+    end
+
+    % Detect Bursts
+    if length(CMAinfo.BurstThreshold) > 1
+        [BurstNumber,BurstStarts,BurstEnds,BurstDurations,NumSpikesInBursts] = ...
+            DetectBurstNoRound(Spikes, CMAinfo.BurstThreshold(s), CMAinfo.TailThreshold(s));
+    else
+        [BurstNumber,BurstStarts,BurstEnds,BurstDurations,NumSpikesInBursts] = ...
+            DetectBurstNoRound(Spikes, CMAinfo.BurstThreshold, CMAinfo.TailThreshold);
+    end
+
+    SpikesIncludedBursts = Spikes(BurstNumber > 0) ; %Index of the all spikes included to bursts
+%     SpikesExcludedBursts = Spikes(BurstNumber == 0) ; %Index of the all spikes excluded from bursts
+
+    % ISI in burst
+    ISIinBurst = BurstDurations ./ (NumSpikesInBursts -1); 
+
+    % Calculate & Collect statistics
+    BI.ID = IDs(s);
+    BI.SpikeCount = length(Spikes);
+    BI.SpikeRate = BI.SpikeCount / simtime;
+    BI.BurstCount = length(BurstStarts);
+    BI.BurstRate = BI.BurstCount / simtime;
+    BI.BurstDuration = mean(BurstDurations);
+    BI.SpikesInBurst = mean(NumSpikesInBursts);
+    BI.BurstSpikeRatio = sum(NumSpikesInBursts) / BI.SpikeCount; %%% can be 0 / nan
+    BI.BurstSpikes = {SpikesIncludedBursts};
+    BI.BurstStarts = {BurstStarts};
+    BI.BurstEnds = {BurstEnds};
+    BI.ISIinBurst = mean(ISIinBurst); 
+    BI.TypeOfSpike = {BurstNumber};
+
+    if s == 1
+        BurstInfo = BI;
+    else
+        BurstInfo = [BurstInfo;BI];
+    end
+
+end
+BurstInfo = struct2table(BurstInfo);
+end
+

DetectBurstNoRound.m

@@ -0,0 +1,119 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% DetectBurstNoRound.m: Detects the bursts using CMA-algorithm
+% authors: Inkeri A. Välkki, Kerstin Lenk, Jarno E. Mikkonen, Fikret E. Kapucu, Jari A. K. Hyttinen
+% date: 2016 - 2019
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [BurstNumbers, BurstStarts, BurstEnds, BurstDurations, NumSpikesInBursts] = DetectBurstNoRound(SpikeTimes, BurstThreshold, TailThreshold)
+% [BurstNumbers, BurstStarts, BurstEnds, BurstDurations, NumSpikesInBursts] = DetectBurstNoRound(SpikeTimes, BurstThreshold, TailThreshold)
+% in:
+%   SpikeTimes: spiketimes in ms, sorted, vector
+%   BurstThreshold, TailThreshold: thresholds for burst detection, scalars
+% out:
+%   BurstNumbers: burst number to which the spike belongs to, 0 for
+%   non-burst spikes
+%   BurstStarts: start times of bursts (ms)
+%   BurstEnds: end times of bursts (ms)
+%   BurstDurations: burst durations (ms)
+%   NumSpikesInBursts: number of spikes in each burst
+
+% Number of spikes
+nspikes = length(SpikeTimes);
+
+% For keeping track of spike types
+TypeOfSpike = zeros(1,nspikes); % 1 for bursts, 2 for tails , 0 for individuals
+
+% Determine burst spikes
+isis = find(diff(SpikeTimes) < BurstThreshold);
+bspikes = false(size(SpikeTimes));
+bspikes(isis) = 1;
+bspikes(isis+1) = 1;
+TypeOfSpike(bspikes) = 1;
+
+%% omitting bursts having fewer spikes then desired (in this case 2 spikes)
+
+%+++comment out if two spike bursts are wanted
+minNOspikes = 3; % specify the min number of spikes in a burst to be
+BreakPointAfter=[];
+BreakPointBefore=[];
+burstspikes = find(TypeOfSpike==1);
+if ~isempty(burstspikes)
+    BreakPoints = find(diff(SpikeTimes(burstspikes)) > BurstThreshold);
+    BreakPointBefore = [burstspikes(1) burstspikes(BreakPoints+1)];
+    BreakPointAfter = [burstspikes(BreakPoints) burstspikes(end)];
+end
+
+ZeroIndexes = find(BreakPointAfter-BreakPointBefore+ 1<minNOspikes); % spike indexes to be nullified
+for l=1:length(ZeroIndexes)
+    TypeOfSpike(BreakPointBefore(ZeroIndexes(l)):BreakPointAfter(ZeroIndexes(l))) = 0;
+end
+%+++make comment out till here++++
+
+%% TAIL CALCULATION
+
+% Find tail spikes
+isis = find(diff(SpikeTimes) < TailThreshold);
+tspikes = false(size(SpikeTimes));
+tspikes(isis) = 1;
+tspikes(isis+1) = 1;
+tspikes(bspikes) = 0;
+TypeOfSpike(tspikes) = 2;
+
+%% omitting tails away from bursts and merging rest with the bursts
+
+tailspikes = find(TypeOfSpike==2);
+controlForLoop = TypeOfSpike;
+resultOfmerging=zeros(1,nspikes);
+%/loop till all the burst related tails merged and counted as bursts
+while ~isequal(controlForLoop,resultOfmerging)
+    controlForLoop = TypeOfSpike;
+    for t = tailspikes
+        if t== 1
+            if (SpikeTimes(t + 1) - SpikeTimes(t) <= TailThreshold...
+                    && TypeOfSpike(t + 1) == 1)
+                TypeOfSpike(t) = 1;
+            end
+        elseif t == nspikes
+            if(SpikeTimes(t) - SpikeTimes(t-1) <= TailThreshold...
+                    && TypeOfSpike(t - 1) == 1)
+                TypeOfSpike(t) = 1;
+            end
+        else
+            if (TypeOfSpike(t + 1) == 1 && SpikeTimes(t + 1) - SpikeTimes(t) <= TailThreshold) ...
+                    || (TypeOfSpike(t - 1) == 1 && SpikeTimes(t) - SpikeTimes(t-1) <= TailThreshold)
+                TypeOfSpike(t) = 1;
+            end
+        end
+    end
+    resultOfmerging = TypeOfSpike;
+    tailspikes = find(TypeOfSpike==2);
+end
+
+nonbursts = TypeOfSpike ~= 1;
+TypeOfSpike(nonbursts) = 0;
+
+%%Burst Starts  and Burst Ends
+burstspikes = find(TypeOfSpike==1);
+if ~isempty(burstspikes)
+    BreakPoints = find(diff(SpikeTimes(burstspikes)) > TailThreshold);
+    BreakPointStarts = [burstspikes(1) burstspikes(BreakPoints+1)];
+    BreakPointEnds = [burstspikes(BreakPoints) burstspikes(end)];
+
+    BurstStarts = SpikeTimes(BreakPointStarts);
+    BurstEnds = SpikeTimes(BreakPointEnds);
+
+    NumSpikesInBursts = BreakPointEnds - BreakPointStarts+1;
+    BurstDurations = BurstEnds - BurstStarts;
+    %     IBI = BurstStarts - BurstEnds;
+    %     ISIinBurst
+    BurstNumbers = zeros(size(TypeOfSpike));
+    for b = 1:length(BreakPointStarts)
+        BurstNumbers(BreakPointStarts(b):BreakPointEnds(b)) = b;
+    end
+else
+    BurstStarts = [];
+    BurstEnds= [];
+    BurstDurations= [];
+    NumSpikesInBursts= [];
+    BurstNumbers = TypeOfSpike;
+end

alphas2ths.m

@@ -0,0 +1,32 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% alphas2ths.m: Calculates the thresholds for burst detection
+% authors: Inkeri A. Välkki, Kerstin Lenk, Jarno E. Mikkonen, Fikret E. Kapucu, Jari A. K. Hyttinen
+% date: 2016 - 2019
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [BurstThreshold, TailThreshold] = alphas2ths(CMAcurve, BurstAlpha, TailAlpha)
+% [BurstThreshold, TailThreshold] = alphas2ths(CMAcurve, BurstAlpha, TailAlpha)
+% in:
+%   CMAcurve: of ISIinfo, vector
+%   BurstAlpha, TailAlpha: scalars
+% out:
+%   BurstThreshold, TailThreshold: scalars
+
+MaxCMA =  max(CMAcurve);
+MaxPointCMA = find(CMAcurve==max(CMAcurve),1, 'last');
+
+if isnan(BurstAlpha) || isempty(CMAcurve)
+    BurstThreshold = nan;
+else
+    BurstDistants = abs(CMAcurve(MaxPointCMA:end) - (BurstAlpha * max(MaxCMA)));
+    BurstThreshold = MaxPointCMA + find(BurstDistants==min(BurstDistants), 1, 'last') -1;
+end
+
+if isnan(TailAlpha) || isempty(CMAcurve)
+    TailThreshold = nan;
+else
+    TailDistants = abs(CMAcurve(BurstThreshold:end) - (TailAlpha * max(MaxCMA)));
+    TailThreshold = BurstThreshold + find(TailDistants==min(TailDistants), 1, 'last') -1;
+end
+
+end

getCMAinfo.m

@@ -0,0 +1,83 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% getCMAinfo.m: Caculates the parameters for CMA burst detection
+% authors: Inkeri A. Välkki, Kerstin Lenk, Jarno E. Mikkonen, Fikret E. Kapucu, Jari A. K. Hyttinen
+% date: 2016 - 2019
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [CMAinfo, ISIinfo] = getCMAinfo(CMAType, SpikeTimes, IDs)
+% [CMAinfo] = getCMAinfo(CMAType, ISIinfo)
+% in:
+%   CMAtype: 'original', 'net' or 'net_alpha'
+%   SpikeTimes: SpikeTimes in ms in a cell
+%   IDs: IDs for the spike sequences
+% out:
+%   CMAinfo: CMA parameters
+%   ISIinfo: ISI statistics
+
+CMAinfo.type = CMAType;
+CMAinfo.ID = IDs;
+
+% Emre's original CMA
+if strcmp(CMAType, 'original')
+
+    % Calculate ISIs & ISIinfo for all the SpikeTimes
+    ISIs = cellfun(@diff, SpikeTimes, 'UniformOutput', false);
+    ISIinfo = getISIinfo(ISIs, IDs);
+
+    % Calculate the alphas
+    [CMAinfo.BurstAlpha, CMAinfo.TailAlpha] = skw2alpha(ISIinfo.skewISI);
+
+    CMAinfo.BurstThreshold = nan(size(SpikeTimes));
+    CMAinfo.TailThreshold = nan(size(SpikeTimes));
+    % Calculate CMA thresholds
+    for n = 1:length(SpikeTimes)
+        [CMAinfo.BurstThreshold(n), CMAinfo.TailThreshold(n)] = alphas2ths(ISIinfo.CMAcurve{n}, CMAinfo.BurstAlpha(n), CMAinfo.TailAlpha(n));
+    end
+
+% Common alpha value for all the neurons !! Is not such a good idea to use
+elseif strcmp(CMAType, 'net_alpha')
+
+    % Calculate ISIs and ISIinfo
+    ISIs = cellfun(@diff, SpikeTimes, 'UniformOutput', false);
+    ISIinfo = getISIinfo(ISIs, IDs);
+
+    % Combined ISI histogram
+    ISIdataAll = horzcat(ISIs{:});
+    ISIinfoAll = getISIinfo({ISIdataAll}, {'all'});
+    ISIinfo = [ISIinfo;ISIinfoAll];
+
+    % Calculate the alphas
+    [CMAinfo.BurstAlpha, CMAinfo.TailAlpha] = skw2alpha(ISIinfoAll.skewISI);
+
+    CMAinfo.BurstThreshold = nan(size(SpikeTimes));
+    CMAinfo.TailThreshold = nan(size(SpikeTimes));
+    % Calculate CMA thresholds
+    for n = 1:length(SpikeTimes)
+        [CMAinfo.BurstThreshold(n), CMAinfo.TailThreshold(n)] = alphas2ths(ISIinfo.CMAcurve{n}, CMAinfo.BurstAlpha, CMAinfo.TailAlpha);
+    end
+
+% Common thresholds for burst detection
+elseif strcmp(CMAType, 'net')
+
+    % Calculate ISIs and ISIinfo
+    ISIs = cellfun(@diff, SpikeTimes, 'UniformOutput', false);
+
+    % Combined ISI histogram
+    ISIdataAll = horzcat(ISIs{:});
+    ISIinfoAll = getISIinfo({ISIdataAll}, {'all'});
+    ISIinfo = ISIinfoAll;
+
+    % Calculate the burst alphas
+    [CMAinfo.BurstAlpha, CMAinfo.TailAlpha] = skw2alpha(ISIinfoAll.skewISI);
+
+    % Calculate CMA thresholds
+    [CMAinfo.BurstThreshold, CMAinfo.TailThreshold] = alphas2ths(ISIinfoAll.CMAcurve{1}, CMAinfo.BurstAlpha, CMAinfo.TailAlpha);
+
+
+end
+
+
+end
+
+
+

getISIinfo.m

@@ -0,0 +1,60 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% getISIinfo.m: Calculates the ISI statistics
+% authors: Inkeri A. Välkki, Kerstin Lenk, Jarno E. Mikkonen, Fikret E. Kapucu, Jari A. K. Hyttinen
+% date: 2016 - 2019
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ISIinfo] = getISIinfo(ISIs, IDs)
+% getISIinfo calculates ISI statistics
+% function [ISIinfo] = getISIinfo(ISIs, IDs)
+% in:
+%   ISIs: cell of ISI sequences in ms
+%   ID: IDs for the ISI sequences
+% out:
+%   ISIinfo: table of ISI statistics
+
+% Calculate ISI statistics for all ISIs
+ISIinfo = cellfun(@calc_stats, ISIs, IDs);
+
+% Convert to table
+ISIinfo = struct2table(ISIinfo);
+
+end
+
+function [info] = calc_stats(ISI, ID)
+% calc_stats calculates II statistics of one ISI squence
+info.ID = ID;
+info.ISIs = ISI;
+
+if isempty(ISI)
+    info.maxISI = nan;
+    info.minISI = nan;
+    info.meanISI = nan;
+    info.stdISI = nan;
+    info.medianISI = nan;
+    info.histISI = [];
+    info.cumhistISI = [];
+    info.CMAcurve = [];
+    info.skewCMA = nan;
+    info.skewISI = nan;
+else
+    info.maxISI = max(ISI);
+    info.minISI = min(ISI);
+    info.meanISI = mean(ISI);
+    info.stdISI = std(ISI);
+    info.medianISI = median(ISI);
+    info.histISI = histc(ISI , (0:1:20000));
+    info.cumhistISI = cumsum(info.histISI);
+    info.CMAcurve = info.cumhistISI ./ (1:1:20001);
+    %SkewnessCMA
+    if  info.maxISI > 20000 || isnan(info.maxISI)
+        info.skewCMA =  skewness(info.CMAcurve);
+    else
+        info.skewCMA = skewness(info.CMAcurve(1:ceil(info.maxISI)));
+    end
+    info.skewISI = skewness(ISI);
+end
+    info.histISI = {info.histISI};
+    info.cumhistISI = {info.cumhistISI};
+    info.CMAcurve = {info.CMAcurve};
+end