function [Bursts, Pauses] = RGSDetect(Spikes, N_min, Steps, p, alpha) |
%% RGSDETECT help |
% |
% Determines burst and pause interspike interval (ISI) thresholds and |
% identifies burst and pause strings based on the Robust Gaussian Surprise |
% (RGS) method. |
% |
% INPUTS: |
% Spikes (s): Times of spikes in seconds. |
% N_min: Minimum number of spikes to be considered a burst/pause |
% string. |
% Steps (log10(s)): Bin edges for histogram count (histc) of the |
% log ISIs. |
% p: Bottom and top p% used as outliers to calculate central |
% location; keep p in range [0.05, 0.30] (default 0.05). |
% alpha: Value used in Bonferroni correction; lower value of |
% alpha to filter out false positives (default 0.05). |
% NOTE: Requires MATLAB statistics toolbox. |
% |
% OUTPUTS: |
% Bursts: Structure containing burst information. |
% Bursts.BurstingSpikes (s): Column of times of all spike times |
% included in a burst. |
% Bursts.IBF (Hz): Column of intraburst frequency (IBF) of each |
% burst. |
% Bursts.NumSpikes: Column of number of spikes in each burst. |
% Bursts.Windows (s): 2 Columns of start and end times of each |
% burst. |
% Pauses: Structure containing pause information. |
% Pauses.AllSpikes (s): Start times of all ISIs that satisfy |
% pause threshold. |
% Pauses.AllLengths (s): Lengths of all ISIs that satisfy pause |
% threshold. |
% Pauses.PausingSpikes (s): Column of all spike times |
% included in a pause string. |
% Pauses.IPF (Hz): Column of intrapause frequency (IPF) of each |
% pause string. |
% Pauses.NumSpikes: Column of number of spikes in each pause |
% string. |
% Pauses.Windows (s): 2 Columns of start and end times of each |
% pause string. |
% NOTE: Rows of structure elements correspond to the same burst or |
% pause. |
% NOTE: Normalized Log ISI Distribution (NLISI) plot is used confirm |
% central distribution is centered on 0. If distribution is not |
% centered on 0, change p until it is. Use steps to adjust the |
% x-axis. |
% |
% EXAMPLE values used in this paper: |
% Spikes = (Spike times go here); |
% N_min = 2; |
% Steps = -3:0.005:1.5; |
% p = 0.05; |
% alpha = 0.05; |
% [Bursts, Pauses] = RGSDetect(Spikes, N_min, Steps, p, alpha); |
% |
% REFERENCE: Ko et al. (2012) |
% Detection of bursts and pauses in spike trains |
% J Neurosci Methods 211:145–158 |
% |
%%% NORMALIZED LOG ISI DISTRIBUTION |
ISIs = Spikes(2:end)-Spikes(1:end-1); %Offset spikes by 1 and subtract |
%for ISI's |
LogISIs = log10(ISIs); %Take the log10 of ISI's |
N = length(LogISIs); %N = number of ISIs |
Q = max([20,floor(0.2×N)]); %Set window length as max of 20 and 20% of N |
NLISITrain = zeros(1,length(LogISIs))'; %Initialize Normalized Log ISI Train |
%Central Location of first 2×Q+1 ISIs |
CentralLoc1 = ComputeCL(LogISIs(1:2×Q+1), Steps, p); |
%Subtract Central Location (Normalize) |
NLISITrain(1:Q) = LogISIs(1:Q) - CentralLoc1; |
%Central Location of last 2×Q+1 ISIs |
CentralLoc2 = ComputeCL(LogISIs(N-2×Q:N), Steps, p); |
%Subtract Central Location (Normalize) |
NLISITrain(N-Q+1:end) = LogISIs(N-Q+1:end) - CentralLoc2; |
%For the middle portion |
for i = Q+1:N-Q |
%Compute central location for portion of index +/− Q and subtract |
NLISITrain(i) = LogISIs(i) - ComputeCL(LogISIs(i-Q:i+Q), Steps, p); |
end |
%Get statistics of the NLISI train |
med = median(NLISITrain); |
pool_MAD = mad(NLISITrain); |
CentralDistBounds = [med - pool_MAD×2.58 med + pool_MAD×2.58]; |
mu = median(NLISITrain); |
sigma = mad(NLISITrain); |
%Plot the NLISI Distribution |
figure |
hold on |
%Run a smoothing pdf kernel. |
NLISIpdf = pdf(fitdist(NLISITrain,'Kernel'), Steps); |
NLISIpdf = NLISIpdf./sum(NLISIpdf); |
plot(Steps, NLISIpdf,'g') |
%Plot treshold lines |
plot([CentralDistBounds(1) CentralDistBounds(1)], [0 max(NLISIpdf)], '–r') |
plot([CentralDistBounds(2) CentralDistBounds(2)], [0 max(NLISIpdf)], '–b') |
xlabel 'Normalized Log ISIs' |
ylabel 'Probability' |
title 'Normalized Log ISI Distribution' |
%% BURST AND PAUSE STRING DETECTION |
%Get index and ISI lengths of all ISIs that satisfy burst threshold |
Burst_Thresh = CentralDistBounds(1); |
BurstINDXS = 1:length(NLISITrain); |
BurstINDXS(NLISITrain >= Burst_Thresh) = []; %Delete all indexes greater than |
%the burst threshold |
if ∼isempty(BurstINDXS) |
%Matrix of all potential burst ISIs and their indexes |
BurstsM = [NLISITrain(NLISITrain < Burst_Thresh)';BurstINDXS]; |
[∼,c] = size(BurstsM); |
Burst_Seed = mat2cell(BurstsM,2,ones(1,c,1)); |
else |
Burst_Seed = {}; |
Bursts = Burst_Seed; |
end |
%Loop through each potential burst ISI (Burst Seed) |
for i = 1:length(Burst_Seed); |
%Go forward and backward from the current burst until both conditions |
%are unsatisfied |
forward = 1; |
backward = 1; |
while forward ∥ backward |
currBurst = cell2mat(Burst_Seed(i)); |
%Go forward 1 ISI |
if forward |
%Set current ISI as end of the current burst |
currSpike = currBurst(:,end); |
if currSpike(2) ∼= length(NLISITrain) |
%q is number of spikes |
[∼,q] = size(currBurst); |
%P1 is probability burst will occur assuming Gaussian |
%distribution with mean, mu×q, and std, sqrt(q)×sigma |
P1 = normcdf(sum(currBurst(1,:)), mu×q, sqrt(q).×sigma); |
testBurst = [currBurst [NLISITrain(currSpike(2)+1);… |
currSpike(2)+1]]; |
%P2 is the same probability with the next ISI added to the |
%burst |
P2 = normcdf(sum(testBurst(1,:)), mu×(q+1), … |
sqrt(q+1).×sigma); |
%If the next ISI increased the probability of the burst |
%occurring |
if P2 >= P1 |
%Stop going forward |
forward = 0; |
else
|
%Otherwise, set the current burst seed to the tested |
%burst |
Burst_Seed{i} = testBurst; |
end |
else |
%Stop going forward if at the end of the ISI train |
forward = 0; |
end |
end |
currBurst = cell2mat(Burst_Seed(i)); |
%Go backward 1 ISI |
if backward |
%Set current ISI as end of the current burst |
currSpike = currBurst(:,1); |
if currSpike(2) ∼= 1 |
%q is number of spikes |
[∼,q] = size(currBurst); |
%P1 is probability burst will occur assuming Gaussian |
%distribution with mean, mu×q, and std, sqrt(q)×sigma |
P1 = normcdf(sum(currBurst(1,:)), mu×q, sqrt(q).×sigma); |
testBurst = [[NLISITrain(currSpike(2)-1);currSpike(2)-1] … |
currBurst]; |
%P2 is the same probability with the next ISI added to the |
%burst |
P2 = normcdf(sum(testBurst(1,:)), mu*(q+1), … |
sqrt(q+1).*sigma); |
%If the next ISI increased the probability of the burst |
%occurring |
if P2 >= P1 |
%Stop going backward |
backward = 0; |
else |
%Otherwise, set the current burst seed to the tested |
%burst |
Burst_Seed{i} = testBurst; |
end |
else |
%Stop going backward if at the end of the ISI train |
backward = 0; |
end |
end |
end |
end |
if ∼isempty(Burst_Seed) |
%Initialize BurstInfo |
BurstInfo = zeros(length(Burst_Seed),3); |
%Get start index of each burst
|
BurstInfo(:,1) = cellfun(@(x) x(2,1),Burst_Seed); |
%Get end index of each burst |
BurstInfo(:,2) = cellfun(@(x) x(2,end),Burst_Seed); |
%Get P-value of each burst (probability of occurence assuming Gaussian |
%distribution) |
BurstInfo(:,3) = cellfun(@(x) normcdf(sum(x(1,:)), mu*length(x), … |
sqrt(length(x)).×sigma),Burst_Seed); |
%Filter out bursts less than minimum number of spikes specified by N_min |
BurstInfo(BurstInfo(:,2)-BurstInfo(:,1)+2 < N_min,:) = []; |
%Filter out overlapping bursts |
no_overlap = 0; |
i=1; |
if ∼isempty(BurstInfo) |
[r,∼] = size(BurstInfo); |
if r ∼= 1 |
while ∼no_overlap |
%If the indexes of the burst ISIs don't intersect |
if isempty(intersect(BurstInfo(i,1):BurstInfo(i,2),… |
BurstInfo(i+1,1):BurstInfo(i+1,2))) |
%move to the next burst |
i = i+1; |
else |
%If the intersect, choose the burst with the lower P |
%value |
if BurstInfo(i,3) <= BurstInfo(i+1,3) |
BurstInfo(i+1,:) = []; |
else |
BurstInfo(i,:) = []; |
end
|
end |
%When the end is reached, stop |
[r,∼] = size(BurstInfo); |
if i == r |
no_overlap = 1; |
end |
end |
end |
%r is the number of rows or the number of bursts |
[r,∼] = size(BurstInfo); |
Bursts.BurstingSpikes = []; |
%for each burst, append the burst spikes |
for i = 1:r |
Bursts.BurstingSpikes = [Bursts.BurstingSpikes;… |
Spikes(BurstInfo(i,1):BurstInfo(i,2)+1)]; |
end |
end |
%Bonferroni correction for false positives |
KB = length(find(BurstInfo(:,3) < alpha)); |
BurstInfo(BurstInfo(:,3)×KB >= alpha,:) = []; |
%Use the indexes in burst info to find the burst windows |
Bursts.Windows = [Spikes(BurstInfo(:,1)) Spikes(BurstInfo(:,2)+1)]; |
%Use the indexes to find the number of spikes in each burst |
Bursts.NumSpikes = BurstInfo(:,2) - BurstInfo(:,1) + 2; |
%Use the number of spikes and windows to calculate the IBF |
Bursts.IBF = Bursts.NumSpikes./(Bursts.Windows(:,2) - … |
Bursts.Windows(:,1)); |
end |
%Get index and ISI lengths of all NLISIs that satisfy pause threshold
|
Pause_Thresh = CentralDistBounds(2); |
PauseINDXS = 1:length(NLISITrain); |
%Delete all indexes less than the pause threshold |
PauseINDXS(NLISITrain <= Pause_Thresh) = []; |
if ∼isempty(PauseINDXS) |
%Matrix of all potential pause string NLISIs and their indexes |
PausesM = [NLISITrain(NLISITrain > Pause_Thresh)';PauseINDXS]; |
[∼,c] = size(PausesM); |
Pause_Seed = mat2cell(PausesM,2,ones(1,c,1)); |
else
|
Pause_Seed = {}; |
Pauses = []; |
end |
%Loop through each potential pause string NLISI (Pause Seed) |
for i = 1:length(Pause_Seed); |
%Go forward and backward from the current pause string until both conditions |
%are unsatisfied |
forward = 1; |
backward = 1; |
while forward ∥ backward |
currPause = cell2mat(Pause_Seed(i)); |
%Go forward 1 ISI |
if forward |
%Set current ISI as end of the current pause string |
currPauseind = currPause(:,end); |
if currPauseind(2) ∼= length(NLISITrain) |
[∼,q] = size(currPause); |
%P1 is probability pause string will occur assuming Gaussian |
%distribution with mean, mu×q, and std, sqrt(q)×sigma |
P1 = (1-normcdf(sum(currPause(1,:)), mu×q, sqrt(q).×sigma));
|
testPause = [currPause [NLISITrain(currPauseind(2)+1);… |
currPauseind(2)+1]]; |
%P2 is the same probability with the next ISI added to the |
%pause string |
P2 = (1-normcdf(sum(testPause(1,:)), mu×(q+1), … |
sqrt(q+1).×sigma)); |
%If the next ISI increased the probability of the pause |
%string occurring |
if P2 >= P1 |
%Stop going forward |
forward = 0; |
else |
%Otherwise, set the current pause seed to the tested |
%pause string |
Pause_Seed{i} = testPause; |
end |
else |
forward = 0; |
end |
end |
currPause = cell2mat(Pause_Seed(i)); |
%Go backward 1 ISI |
if backward |
currPauseind = currPause(:,1); |
if currPauseind(2) ∼= 1 |
[∼,q] = size(currPause); |
%P1 is probability pause string will occur assuming Gaussian |
%distribution with mean, mu×q, and std, sqrt(q)×sigma |
P1 = (1-normcdf(sum(currPause(1,:)), mu×q, sqrt(q).×sigma)); |
testPause = [[NLISITrain(currPauseind(2)-1);… |
currPauseind(2)-1] currPause]; |
%P2 is the same probability with the next ISI added to the |
%pause string |
P2 = (1-normcdf(sum(currPause(1,:)), mu×(q+1), … |
sqrt(q+1).×sigma)); |
%If the next ISI increased the probability of the pause |
%string occurring |
if P2 >= P1 |
%Stop going forward |
backward = 0; |
else |
%Otherwise, set the current pause seed to the tested |
%pause string |
Pause_Seed{i} = testPause; |
end |
else |
backward = 0; |
end |
end |
end |
end |
if ∼isempty(Pause_Seed) |
%Initialize PauseInfo variable |
PauseInfo = zeros(length(Pause_Seed),3); |
%Starting indexes of pause strings |
PauseInfo(:,1) = cellfun(@(x) x(2,1),Pause_Seed); |
%Ending indexes of pause strings |
PauseInfo(:,2) = cellfun(@(x) x(2,end),Pause_Seed); |
%P-value of the pause strings |
PauseInfo(:,3) = cellfun(@(x) (1-normcdf(sum(x(1,:)), mu×length(x), … |
sqrt(length(x)).×sigma)),Pause_Seed); |
%Minimum number of spikes filter |
PauseInfo(PauseInfo(:,2)-PauseInfo(:,1) + 2 < N_min,:) = []; |
%Filter out overlaps |
no_overlap = 0; |
i=1; |
if ∼isempty(PauseInfo) |
%r is number of current pause strings |
[r,∼] = size(PauseInfo); |
if r ∼= 1 |
while ∼no_overlap |
%If the indexes of the burst ISIs don't intersect |
if isempty(intersect(PauseInfo(i,1):PauseInfo(i,2),… |
PauseInfo(i+1,1):PauseInfo(i+1,2))) |
%Move to next pause string |
i = i+1; |
else |
%Choose the pause string with lower P-value |
if PauseInfo(i,3) <= PauseInfo(i+1,3) |
PauseInfo(i+1,:) = []; |
else |
PauseInfo(i,:) = []; |
end |
end |
%End if the last pause string is reached |
[r,∼] = size(PauseInfo); |
if i == r |
no_overlap = 1; |
end |
end |
end |
end |
%Use indexes to find start and end times |
Pauses.Windows = [Spikes(PauseInfo(:,1)) Spikes(PauseInfo(:,2)+1)]; |
%Use indexes to determine number of spikes |
Pauses.NumSpikes = PauseInfo(:,2) - PauseInfo(:,1) + 2; |
%Use windows and numspikes to calculate IPF |
Pauses.IPF = Pauses.NumSpikes./(Pauses.Windows(:,2) - … |
Pauses.Windows(:,1)); |
Pauses.PausingSpikes = []; |
%Add pausing spikes from each pause string to the pausingspikes element |
[r,∼] = size(PauseInfo); |
for i = 1:r |
Pauses.PausingSpikes = [Pauses.PausingSpikes;… |
Spikes(PauseInfo(i,1):PauseInfo(i,2)+1)]; |
end |
end |
%Bonferroni correction |
KP = length(find(PauseInfo(:,3) < alpha)); |
PauseInfo(PauseInfo(:,3)×KP >= alpha,:) = []; |
%Get all pauses using the pause indexes |
Pauses.AllSpikes = Spikes(PauseINDXS); |
%Get the lengths of all the pauses that satisfy the threshold |
Pauses.AllLengths = ISIs(PauseINDXS); |
end
|