function eRT(testtrace1,prewin11,dt1, templ_pos_from_boundary, increment1)
% This function calculates the effective recovery time (eRT) and the
% effective synaptic integration time (eSIT) for an input trace
% "testtrace1" which has spikes in it. The testtrace1 is e.g. current
% clamp data of intracellular recording from a neuron.
%
% Prewin1 is the size of the window symmetrically around the spike, in
% miliseconds. Increment is the distance in ms from each testpoint.
% dt1 is the sampling interval, templ_pos_from_boundary is
% the position of test template from boundary.
% For more information please see Berg et al 2008 PLoS One:
% http://dx.plos.org/10.1371/journal.pone.0003218
% Suggested values for parameters:
%
% prewin11= 40 ms
% templ_pos_from_boundary = 3 ms
% increment = 0.5 ms
increment = increment1/(dt1*1000);
prewin1 = prewin11/(dt1*1000);
startpoint= 1;
prewin=startpoint:increment:prewin1;
postwin=prewin;
clear eRT_A; clear eSIT_A; clear num_spikes_A
for ii=1:length(prewin)
% Getting the effecttive Recovery time
spikeava=0;spikeavb=0;
[outs1b outs2b]=spike_select(testtrace1,prewin(ii),postwin(ii));
if outs1b(1) <1, break, end
spikeavb=spikeavb+outs1b; nspikesb1=spikeavb(1);nspikestotalb=nspikesb1;
timebasespikeb= -prewin(ii)*dt1:dt1:(size(outs2b,2)-prewin(ii)-1)*dt1;
loc_pre=max(-0.01,-prewin(ii)*dt1);loc_post=min(0.03,postwin(ii)*dt1);
loc_pre_cord=min(find(timebasespikeb > loc_pre));loc_post_cord=max(find(timebasespikeb < loc_post));
if size(outs2b,2) < 900
loc_post_cord=size(outs2b,2);
else
loc_post_cord=900;
end
templatedist1=outs2b(:,:)';
limit_template=200;
jj=0;
loc_pre_cord=jj+templ_pos_from_boundary/(1000*dt1); %adding templ_pos_from_boundary to stay clear of previous spikes.
if size(templatedist1,1) > loc_pre_cord
templatedist=reshape(templatedist1(loc_pre_cord,:),1,[]);
else
templatedist=reshape(templatedist1(jj+1,:),1,[]);
end
clear h_a; clear p_a
for kk=1:size(outs2b,2)
testdist=outs2b(:,kk)';
[h_post,p_post,k_post] = kstest2(reshape(testdist,1,[]),templatedist);
h_a(kk)=h_post;
p_a(kk)=p_post;
h_locus(kk)=kk;
tball_test(kk)=kk;
end
%location for the effective RT : Choosing outside +- 0.5 ms
av_trace=mean(outs2b',2);
Resetpotential= find(av_trace(find(timebasespikeb>0))== min(av_trace(find(timebasespikeb>0))))...
+min(find(timebasespikeb>0))-1;
loc_RT_range= Resetpotential:size(outs2b,2) ;
locybase=ones(1,length(loc_RT_range));
%plot(loc_RT_range,locybase,'.')
% loc_RT_range=find( timebasespikeb> 0.001);
loc_SIT_range=find( timebasespikeb < -0.001);
%h_a(loc_RT_range);
% Location of eRT
if isempty(min(find(h_a(loc_RT_range) < 1 ))) < 1 && (length(timebasespikeb) > (min(find(h_a(loc_RT_range) < 1 ))+loc_RT_range(1)))
postspike_change_A = timebasespikeb(min(find(h_a(loc_RT_range) < 1 ))+loc_RT_range(1));
else
postspike_change_A = timebasespikeb(end);
end
% Location of eSIT
if isempty(min(find(h_a(loc_SIT_range) < 1 ))) <1
prespike_change_A = timebasespikeb(max(find(h_a(loc_SIT_range) < 1 )));
else
prespike_change_A =timebasespikeb(1);
end
eRT_A(ii)=postspike_change_A;
eSIT_A(ii)=prespike_change_A;
figure(24)
subplot(1,1,1);plot(timebasespikeb*1000,outs2b','color',[.5 0.5 .5],'LineWidth',.1)
desc = { ['KS test post= ' int2str(h_post)],['p= ' num2str(p_post)]};
text(2,25,['eRT = ' num2str(postspike_change_A*1000) ' ms'])
text(-.03*1000,15,['eSIT ave = ' num2str(prespike_change_A*1000) ' ms'])
text(-.03*1000,0,'KS test ')
text(-.03*1000,-10,'P-value ')
hold on
subplot(1,1,1);plot(timebasespikeb*1000,mean(outs2b',2),'color',[1 0.5 .5],'LineWidth',.3)
plot(timebasespikeb*1000, h_a'*10,'-')
plot(timebasespikeb*1000, p_a'*10-20,'-')
plot(prespike_change_A*1000, 20 ,'o',postspike_change_A, 20 ,'o')
grid on
axis([min(timebasespikeb)*1.2*1000 max(timebasespikeb)*1.2*1000 -80 50])
%axis([0 length(outs2b(1,:)) -75 25])
title([' Number of spikes = ' int2str(outs1b(1)) '/' int2str(outs1b(2))])
xlabel('Time [msec]')
ylabel('Membrane potential [mV]')
hold off
num_spikes_A(ii)=outs1b(1);
end
figure(10)
subplot(2,1,1);plot(prewin(1:length(eRT_A))*dt1*1000,eRT_A*1000,'o',prewin(1:length(eSIT_A))*dt1*1000,eSIT_A*1000,'o')
title('eRT and eSIT as a function of template locatio')
xlabel('Location of template before spike [ms]')
ylabel('eRT and eSIT [ms]')
grid on
subplot(2,1,2);plot(prewin(1:length(num_spikes_A))*dt1*1000,num_spikes_A,'o')
title('Number of spikes in distrubution')
xlabel('Location of template before spike [ms]')
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [out1 out2] = spike_select(testtrace1,prewin1,postwin1)
% This function calculates the average of all the spikes in the input trace
% "testtrace1" in a window "prewin" before the spike and "postwin" after
% the spike. If there should be another spike occuring with that window
% "postwin" following the spike, that spike is not accounted for.
% The first number in the out1 array is the number of spikes in the
% average. The second number is the total number of spikes detected.
% The following arrays is the average spike.
%
% Rune W. Berg 2005
threshold1=0.5; %The threshold for spike selection, amount multiplied by the peak value. Everything below is disregarded.
if isempty(find(testtrace1>0)) >0
out2 = 0;
out1(1)=0;
out1(2)=0;
out1(3:4)=0;
disp(' No spikes in trace')
else
gim=testtrace1-mean(testtrace1);
clear('set_crossgi')
set_crossgi=find(gim(1:end) > threshold1*max(gim)) ; % setting the threshold
clear('index_shift_neggi');clear('index_shift_pos');
clear('set_cross_plusgi');clear('set_cross_minus')
clear('index_shift_posgi');clear('index_shift_neggi')
index_shift_posgi(1)=min(set_crossgi);
index_shift_neggi(length(set_crossgi))=max(set_crossgi);
for i=1:length(set_crossgi)-1
if set_crossgi(i+1) > set_crossgi(i)+1 ;
index_shift_posgi(i+1)=i;
index_shift_neggi(i)=i;
end
end
%These are the coords in the nerve based smooothing:
set_cross_plusgi= set_crossgi(find(index_shift_posgi)); % find(x) returns nonzero arguments.
set_cross_minusgi= set_crossgi(find(index_shift_neggi)); % find(x) returns nonzero arguments.
set_cross_minusgi(length(set_cross_plusgi))= set_crossgi(end);
nspikes= length(set_cross_plusgi); % Number of pulses, i.e. number of windows.
clear('y');clear('y2')
base1gi=0:1:length(gim)-1;
y(1:nspikes)=0;
base2gi=0:1:length(set_crossgi)-1;
y2(1:length(set_crossgi))=0;
base3gi=0:1:length(gim)-1;
y(1:nspikes)=0;
%%% Making sure that the beginning and end are inside the whole window
% avoiding the endpoint problems:
subtrac=nspikes-min(find(set_cross_minusgi > length(testtrace1)-postwin1 ))+1 ;
if isempty(subtrac) > 0
subtrac1=0;
else
subtrac1=subtrac;
end
av_spike=0.;
count1=0;
% For the first spike:
i=1;
% Finding location of spike
spikemax(i) = min(find(gim(set_cross_plusgi(i):set_cross_minusgi(i)) == ...
max(gim(set_cross_plusgi(i):set_cross_minusgi(i))))) +set_cross_plusgi(i)-1;
% Finding location of following spike to avoid overlap
if length(set_cross_minusgi) > 1
spikemax2(i) = min(find(gim(set_cross_plusgi(i+1):set_cross_minusgi(i+1)) == ...
max(gim(set_cross_plusgi(i+1):set_cross_minusgi(i+1))))) +set_cross_plusgi(i+1)-1;
else % in case only one spike in trace
spikemax2(i) = length(gim)-1;
end
% Finding location of previous spike to avoid overlap
spikemax_prior(i)=1;
if (spikemax2(i) - spikemax(i) > postwin1) && (spikemax(i) - spikemax_prior(i) > prewin1)
if spikemax(i) > prewin1 % making sure no spike is too clos to boundary
av_spike=av_spike+testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
indi_spikes(count1+1,1:length(av_spike))=testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
end
count1=count1 +1;
end % END of first spike analysis
for i=2: nspikes - subtrac1-1
% Finding location of spike
spikemax(i) = min(find(gim(set_cross_plusgi(i):set_cross_minusgi(i)) == ...
max(gim(set_cross_plusgi(i):set_cross_minusgi(i))))) +set_cross_plusgi(i)-1;
% Finding location of following spike to avoid overlap
spikemax2(i) = min(find(gim(set_cross_plusgi(i+1):set_cross_minusgi(i+1)) == ...
max(gim(set_cross_plusgi(i+1):set_cross_minusgi(i+1))))) +set_cross_plusgi(i+1)-1;
% Finding location of previous spike to avoid overlap
spikemax_prior(i)=min(find(gim(set_cross_plusgi(i-1):set_cross_minusgi(i-1)) == ...
max(gim(set_cross_plusgi(i-1):set_cross_minusgi(i-1))))) +set_cross_plusgi(i-1)-1;
if (spikemax2(i) - spikemax(i) > postwin1) && (spikemax(i) - spikemax_prior(i) > prewin1)
if spikemax(i) > prewin1 % making sure no spike is too clos to boundary
av_spike=av_spike+testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
indi_spikes(count1+1,1:length(av_spike))=testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
% plot(testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1)))
end
count1=count1 +1;
end
end
% For the LAST spike:
i=(nspikes - subtrac1);
% Finding location of spike
spikemax(i) = min(find(gim(set_cross_plusgi(i):set_cross_minusgi(i)) == ...
max(gim(set_cross_plusgi(i):set_cross_minusgi(i))))) +set_cross_plusgi(i)-1;
% Finding location of following spike to avoid overlap - just
% choosing last point in trace
spikemax2(i) = length(gim)-1;
% Finding location of previous spike to avoid overlap
if length(set_cross_plusgi) >1
spikemax_prior(i)=min(find(gim(set_cross_plusgi(i-1):set_cross_minusgi(i-1)) == ...
max(gim(set_cross_plusgi(i-1):set_cross_minusgi(i-1))))) +set_cross_plusgi(i-1)-1;
else % in case only one spike in trace:
spikemax_prior(i)=1;
end
if (spikemax2(i) - spikemax(i) > postwin1) && (spikemax(i) - spikemax_prior(i) > prewin1)
if spikemax(i) > prewin1 % making sure no spike is too clos to boundary
av_spike=av_spike+testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
indi_spikes(count1+1,1:length(av_spike))=testtrace1(round(spikemax(i) - prewin1):round(spikemax(i) + postwin1));
end
count1=count1 +1;
end % END of last spike processing
av_spiken=av_spike/count1;
out1(1)=count1;
out1(2)=nspikes;
out1(3:length(av_spiken)+2)=av_spiken;
if count1 == 0
disp(' No spikes present in trace !! ')
out2 = 0;
else
out2=indi_spikes;
end
end
