% This function was originally written for Rodriguez, et al., Cytosolic PSD-95
% Interactor Regulates AMPA Receptor Signaling Independent of PSD-95
% Binding. Network Neurosci, 2020.
% This function takes in the spiketimes and the length of the simulation to
% calculate the "network spike" features: the average inter-network-spike
% interval, the standard dev of the INSI, the number of network spikes for
% the simulation, the network spikes per second, and the coefficient of
% variation for the network spikes.
% This function was written by Erin D. Anderson and can be
% accessed at https://www.seas.upenn.edu/~molneuro/
% Last Updated: 11/14/2023
function [insi_avg,insi_std,NScount,NSrate,CV,nSpikesPerBurst_avg,nSpikesPerBurst_std,burstWidthAvg,burstWidthStd] = ...
INSIbyElectrode(spikeindexes,spiketimes,SimTimeInSeconds,StabilizationTime,nUnits)
if nargin == 2
StabilizationTime = 10;
warning('Setting Stabilization Time to 10 s')
end
thr = 4; % need at least 4 events/bin to call it a network spike
timeBin = 100e-3; % [s] - thr and timeBin derived from Rodriguez thesis
insi_avg = zeros(nUnits,1);
insi_std = zeros(nUnits,1);
NScount = zeros(nUnits,1);
NSrate = zeros(nUnits,1);
CV = zeros(nUnits,1);
nSpikesPerBurst_avg = zeros(nUnits,1);
nSpikesPerBurst_std = zeros(nUnits,1);
burstWidthAvg = zeros(nUnits,1);
burstWidthStd = zeros(nUnits,1);
for ii = 1:nUnits
spiketimes_post = spiketimes(spiketimes(spikeindexes == ii) > StabilizationTime); % spiketimes that happen
% after the settling period for this particular neuron
if numel(spiketimes_post) > 1 && spiketimes_post(end)-spiketimes_post(1) > timeBin
% divide the spiketimes into bins to see if reach threshold for network spike
h = histcounts(spiketimes_post,'BinEdges',...
[spiketimes_post(1):timeBin:spiketimes_post(end)]);
rate = h/60.0/timeBin; % h[0] is the counts and h[1] is the bin edges
i_start = find((rate(2:end)>=thr).*(rate(1:end-1)<thr)>0); % these are
% indexes in the
% histcounts
% bins
% find all the indices > threshold, then check the index
% before, make sure it's 0, then multiply to get the start
% index
i_end = find((rate(1:end-1)>=thr).*(rate(2:end)<thr)>0);
if length(i_end) > length(i_start)
i_end = i_end(2:end); % in case it thinks index = 1 is the end of a burst
elseif length(i_start) > length(i_end)
i_end(end+1) = i_start(end); % in case it thinks the last index is the start of a burst
end
if numel(i_start) <= 1
insi = NaN; % if there are no network spikes
burstWidth = NaN;
else
insi = ( i_start(2:end) - i_start(1:end-1) )*timeBin; % convert to time intervals
burstWidth = (i_end-i_start)*timeBin;
end
insi_avg(ii) = mean(insi);
insi_std(ii) = std(insi);
NScount(ii) = length(i_start);
NSrate(ii) = NScount(ii)/SimTimeInSeconds;
CV(ii) = insi_std(ii)/insi_avg(ii);
burstWidthAvg(ii) = mean(burstWidth);
burstWidthStd(ii) = std(burstWidth);
if NScount(ii) ~=0
nSpikesPerBurst_avg(ii) = mean(h(i_start+1)); % i_start is the onset and i_start+1 is peak
nSpikesPerBurst_std(ii) = std(h(i_start+1));
else
nSpikesPerBurst_avg(ii) = 0;
nSpikesPerBurst_std(ii) = 0;
end
else
insi_avg(ii) = 0;
insi_std(ii) = 0;
NScount(ii) = 0;
NSrate(ii) = 0;
CV(ii) = 0;
nSpikesPerBurst_avg(ii) = 0;
nSpikesPerBurst_std(ii) = 0;
burstWidthAvg(ii) = 0;
burstWidthStd(ii) = 0;
end
end
end