Revision 1732ada9c8b7f1432fe3fd494fa00a27e29425ff authored by simkind on 16 March 2024, 19:13:22 UTC, committed by GitHub on 16 March 2024, 19:13:22 UTC
1 parent 7a7d099
PostBurstAHPanalysis.m
% startsearch = peak_times(1,1)/(si/1000);
function [output] = PostBurstAHPanalysis(data,si,bl,startsearch)
numconsec = 50/(si/1000); % 50ms
belowblidx = find(data(startsearch:end) < bl); % find below BL starting from first spike
if ~isempty(belowblidx)
diffbtwn = diff(belowblidx); % if consecutive bins are below threshold, diff should be 1 for each
AHPstart = [];
counter = 1;
if length(diffbtwn) > numconsec
while isempty(AHPstart)
if sum(diffbtwn(counter:counter+numconsec-1)) == numconsec
AHPstart = counter
end
counter = counter + 1;
if counter >= length(diffbtwn)-numconsec
break
end
end
else
AHPstart = [];
end
if ~isempty(AHPstart)
AHPstart_time = startsearch+belowblidx(AHPstart);
AHPstart_amplitude = data(AHPstart_time);
[AHPpeak_amplitude, C] = min(data(AHPstart_time:end));
AHPpeak_time = C+AHPstart_time;
% Find 1sAHP
onesec = 1000/(si/1000);
AHP1s_time = AHPstart_time+onesec;
AHP1s_meanamplitude = mean(data(AHP1s_time-5:AHP1s_time+5));
Output.AHPstart_time = AHPstart_time*(si/1000);
Output.AHPstart_amplitude = AHPstart_amplitude;
Output.AHPpeak_negative_amplitude = AHPpeak_amplitude;
Output.AHPpeak_negative_time = AHPpeak_time*(si/1000);
Output.AHP1s_time = AHP1s_time*(si/1000);
Output.AHP1s_mean_amplitude = AHP1s_meanamplitude;
else
Output.AHPstart_time = [];
Output.AHPstart_amplitude = [];
Output.AHPpeak_negative_amplitude = [];
Output.AHPpeak_negative_time = [];
Output.AHP1s_time = [];
Output.AHP1s_mean_amplitude = [];
end
else
Output.AHPstart_time = [];
Output.AHPstart_amplitude = [];
Output.AHPpeak_negative_amplitude = [];
Output.AHPpeak_negative_time = [];
Output.AHP1s_time = [];
Output.AHP1s_mean_amplitude = [];
end
%% now find when it returns to baseline - slow AHP
if ~isempty(AHPstart_time)
numconsecreturn = 10/(si/1000);
returnAHPbelowblidx = find(data(AHPstart_time:end) >= bl);
if ~isempty(returnAHPbelowblidx)
returnAHPdiffbtwn = diff(returnAHPbelowblidx); % if consecutive bins are below threshold, diff should be 1 for each
returnAHPpossible = [];
counter = 1;
while isempty(returnAHPpossible)
if sum(returnAHPdiffbtwn(counter:counter+numconsecreturn-1)) == numconsecreturn
returnAHPpossible = counter
end
counter = counter + 1;
if counter >= length(returnAHPdiffbtwn)-numconsecreturn
break
end
end
if ~isempty(returnAHPpossible)
AHPend_time = AHPstart_time+returnAHPbelowblidx(returnAHPpossible);
AHPend_amplitude = data(AHPend_time);
AHPduration = AHPend_time - AHPstart_time;
AHP_AUC = trapz(abs(data(AHPstart_time:AHPend_time))); % area under the curve using trapezoidal numerical integration
else
AHPend_time = [];
AHPend_amplitude = [];
AHPduration = [];
AHP_AUC = [];
end
Results.AHPend_time = AHPend_time*(si/1000);
Results.AHPend_amplitude = AHPend_amplitude;
Results.AHPduration = AHPduration*(si/1000);
Results.AHPpeak = Results.baseline_potential - Results.AHPpeak_negative_amplitude
Results.AHPslow = Results.baseline_potential - Results.AHP1s_mean_amplitude
Results.AHP_AreaUnderCurve = (AHP_AUC*(si/1000))/1000;
end
blline = zeros(1,length(data));
blline(:) = bl;
figure,
plot(1:length(data),data,'b',loc,pks,'r*')
hold on
plot(blline)
xx = zeros(1,length(data));
xx(1,AHPstart_time) = -100;
xx(1,AHPpeak_time) = -100;
xx(1,AHPend_time) = -100;
bar(xx)
% xlim([0 AHPend_time+
% xlim(xminmax)
xlim([0 length(data)])
xlabel('points (multiply si/1000 for ms)')
% set(gca,'xticklabel',num2str(get(gca,'xtick')'))
% set(gca,'yticklabel',num2str(get(gca,'ytick')'))
end
%% Find the time constant
timeconstantthresh = AHPpeak_amplitude + Results.AHPpeak*.632;
numconsec = 50;
% while isempty(taupoint)
belowtau = [];
belowtau = find(timeconstantthresh <= data(AHPpeak_time:end))+ AHPpeak_time-1 ;
belowtau(:,2) = data(belowtau(:,1));
if ~isempty(belowtau)
diffbtwntau = diff(belowtau(:,1));
AHPtau = [];
counter = 1;
while isempty(AHPtau)
if sum(diffbtwntau(counter:counter+numconsec-1)) == numconsec
AHPtau = belowtau(counter,1);
end
counter = counter + 1;
if counter >= length(diffbtwntau)-numconsec
break
end
end
end
Results.TimeConstant_Threshold = timeconstantthresh;
Results.TimeConstant_Location = AHPtau*(si/1000);
Results.TimeConstant_Amplitude = data(AHPtau);
Results.TimeConstant = AHPtau - AHPpeak_time;
Computing file changes ...