function varargout = HFblock_th(varargin)
%find the high frequency blocking threshold for the below model file.
%--------------------------------------------------------------------------
%create('axon');
%global axon
%axon.nseg =41; % change the
%number of segments
%axon.dia =10; % um diameter
%axon.Lmode =2;
%axon.Linter_D =100;
%axon.L =axon.Linter_D*axon.dia*(axon.nseg); %um change the length of axon
%axon.da_D =0.6;
%axon.lnodal =1;
%axon.dn_D =0.33;
%axon.xyzi =[1 0 0]; % unitless
%axon.xyzc =[0 0 0]; % um: the coordinate of the node at the center
%insert(axon,'schwarz')
%Iinstim(axon,1,0.005,500,20000) %intracellular stimulation [variable positon(0-1) amp(uA) pw(us) delay(us)]
%Iexstim([ 0 0 200],'Isinth') % [xyz(um),amp(uA),dur(us),delay(us)]
%--------------------------------------------------------------------------
global zeng zeng2 blockth
if nargin==0
action='initial';
else
action=varargin{1};
end
switch action
case 'initial'
%model setup-----------------------------------------------------
[Filename Path]=uigetfile([zeng2.cd zeng2.FILESEP 'model' zeng2.FILESEP '*.m']);
if strcmp('0',num2str(Filename))
return
end
zexst('open model file',Filename);
%/model setup-----------------------------------------------------
zeng2.options.vestim=1; %1=calculated, 2=imported
zexst('change options','vestim',1);
%------------------------------------
if length(zeng2.varlist)>1
zexst('err','Sorry.. This version support only one "axon" in the model file.');
return
else
for tempzeng=1:length(zeng2.varlist)
feval('global',zeng2.varlist{tempzeng})
end
end
%------------------------------------
zeng.dtmode=1;
zeng.dt=5; %us
zeng.tf=30; %ms
zeng2.rundisp =0;
%------------------------------------
Iexstim=zeng.Iexstim;
Iinstim=zeng.Iinstim;
%------------------------------------
blockth.detect_start=round(Iinstim.delay/zeng.dt*1000);
blockth.cal.Ierrpercent = 0.5; %the maximum percent err for the calculated threshold
blockth.Istart = 1000;
blockth.Ilimit = 50*blockth.Istart;
blockth.freq_all = 10000:-1000:4000;
blockth.freq_num = length(blockth.freq_all);
blockth.delay_num = 8;
blockth.delay_all =zeros(blockth.freq_num,blockth.delay_num)
%------------------------------------
if 0
blockth.cal.Dia=10;
else
blockth.cal.Dia=1:1:20;
%blockth.cal.Dia=5:1:15;
%blockth.cal.Dia=[14 15 21:25];
%blockth.cal.Dia=5:5:10;
end
blockth.cal.Dia_length=length(blockth.cal.Dia);
%-----------------------------------------
blockth.cal.X_number=1;
blockth.cal.X=zeros(blockth.cal.X_number,blockth.cal.Dia_length);
if 0 %systematic 0-L
blockth.cal.X(:,1)=linspace(0,1,blockth.cal.X_number)';
for i=2:blockth.cal.Dia_length
blockth.cal.X(:,i)=blockth.cal.X(:,1);
end
elseif 0 %systematic 0-L/2
blockth.cal.X(:,1)=linspace(0,0.5,blockth.cal.X_number)';
for i=2:blockth.cal.Dia_length
blockth.cal.X(:,i)=blockth.cal.X(:,1);
end
elseif 0 %random
%this will not work unless cal.X has 4 dimensions (for fre,delay, ,)
blockth.cal.X=rand(blockth.cal.X_number,blockth.cal.Dia_length)*2-1; %[-1 +1]
end
eval(['Linter_D=' zeng2.var{1}.name '.Linter_D;'])
for i=1:blockth.cal.Dia_length
blockth.cal.X(:,i)=blockth.cal.X(:,i)*Linter_D*blockth.cal.Dia(i);
end
%Don't sort(X), do sort(Z) only.
if 0 % random_z
%this will not work unless cal.Z has 4 dimensions (for fre,delay, ,)
temp=[50 200]; %range (um)
blockth.cal.Z=rand(blockth.cal.X_number,blockth.cal.Dia_length)*diff(temp)+temp(1);
blockth.cal.Z=sort(blockth.cal.Z);
elseif 0
%for systematic Z
else
blockth.cal.Z=[];
end
%-----------------------------------------
blockth.cal.Ith=zeros(blockth.freq_num,blockth.delay_num,blockth.cal.X_number,blockth.cal.Dia_length);
Total_axons=num2str(blockth.cal.Dia_length*blockth.cal.X_number);
%------------------------------------
for freq_i=1:blockth.freq_num;
blockth.freq = blockth.freq_all(freq_i); %this one is used in Isinth
blockth.delay_all(freq_i,:)=Iinstim.delay+[0:(blockth.delay_num-1)]/blockth.delay_num*1000/blockth.freq;
for delay_i=1:blockth.delay_num
zeng.Iinstim.delay=blockth.delay_all(freq_i,delay_i);
disp(['Fre=' num2str(blockth.freq) ' delay=' num2str(zeng.Iinstim.delay)])
%-----------------------------------------
Total_axons_i=0;
%-----------------------------------------
tic
if 1 %for early termination
Early_Ter_check=1;
Early_Ter=figure(99);
text(0.25,0.5,'close me to stop')
else
Early_Ter_check=0;
end
for i=1:blockth.cal.Dia_length
%-----------------------------------------
if 0 %vary L
temp.L=400000; %um
eval([zeng2.var{1}.name '.dia =blockth.cal.Dia(i);'])
eval([zeng2.var{1}.name '.nseg =1+2*round(temp.L/2/blockth.cal.Dia(i)/' zeng2.var{1}.name '.Linter_D)']);
eval([zeng2.var{1}.name '.L =' zeng2.var{1}.name '.nseg*blockth.cal.Dia(i)*' zeng2.var{1}.name '.Linter_D']);
else %donot alter L
eval([zeng2.var{1}.name '.dia =blockth.cal.Dia(i);'])
end
%-----------------------------------------
for j=1:blockth.cal.X_number
Total_axons_i=Total_axons_i+1;
disp([num2str(Total_axons_i) '/' Total_axons])
disp([' Fre=' num2str(blockth.freq) ' delay=' num2str(zeng.Iinstim.delay) ' Dia=' num2str(blockth.cal.Dia(i)) ' x=' num2str(blockth.cal.X(j,i)) ])
%-------------------------------
for k=1:length(zeng.Iexstim)
zeng.Iexstim(k).xyz(1)=Iexstim(k).xyz(1)+blockth.cal.X(j,i);
end
%-------------------------------
Iminmaxmean=blockth.Istart;
Iminmax(1:2)=Iminmaxmean;
zeng2.amp =Iminmaxmean;
zexst('calculate');
if isempty(find(axon.vm(1,[blockth.detect_start:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
%-------------------------------
disp([' start with I = ' num2str([zeng2.amp]) ' - ' num2str(Blocked)])
if Blocked==0
while Blocked==0
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmax(1)=Iminmax(2);
Iminmax(2)=Iminmax(2)*2;
zeng2.amp =Iminmax(2);
zexst('calculate');
if isempty(find(axon.vm(1,[blockth.detect_start:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
disp([' I = ' num2str([zeng2.amp]) ' - ' num2str(Blocked)])
if Iminmax(2)>blockth.Ilimit
Blocked=1
end
end
elseif Blocked==1
while Blocked==1
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmax(2)=Iminmax(1);
Iminmax(1)=Iminmax(1)*0.75;
zeng2.amp =Iminmax(1);
zexst('calculate');
if isempty(find(axon.vm(1,[blockth.detect_start:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
disp([' I = ' num2str([zeng2.amp]) ' - ' num2str(Blocked)])
end
else
disp('err')
return
end
% disp([' I found you at ' num2str(Iminmax)])
while 1
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmaxmean=mean(Iminmax);
if abs(Iminmaxmean)>blockth.Ilimit
blockth.cal.Ith(freq_i,delay_i,j,i)=blockth.Ilimit;
break
end
zeng2.amp =Iminmaxmean;
zexst('calculate');
if isempty(find(axon.vm(1,[blockth.detect_start:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
disp([' I = ' num2str(zeng2.amp) ' - ' num2str(Blocked)])
if Blocked==1
if (100*(Iminmaxmean-Iminmax(1))/Iminmaxmean)>blockth.cal.Ierrpercent
Iminmax(2)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmaxmean])])
blockth.cal.Ith(freq_i,delay_i,j,i)=Iminmaxmean;
break
end
elseif Blocked==0
if (100*(Iminmax(2)-Iminmaxmean)/Iminmaxmean)>blockth.cal.Ierrpercent
Iminmax(1)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmax(2)])])
blockth.cal.Ith(freq_i,delay_i,j,i)=Iminmax(2);
break
end
end
end
end %for j=1:blockth.cal.X_number
HFblock_th('Save Threshold')
end% for i=1:blockth.cal.Dia_length
HFblock_th('Save Threshold')
temp=toc;
if 0
figure
if size(blockth.cal.Ith(freq_i,delay_i,:,:),1)==1
plot(blockth.cal.Dia,blockth.cal.Ith(freq_i,delay_i,:,:));
else
errorbar(blockth.cal.Dia,mean(blockth.cal.Ith(freq_i,delay_i,:,:)),std(blockth.cal.Ith(freq_i,delay_i,:,:)),'b');
end
end
if temp<=60
disp([' finish run with ' num2str(temp) ' seconds' ])
elseif temp<=3600
temp_minute=fix(temp/60);
temp_second=mod(temp,60);
disp([' finish run with ' num2str(temp_minute) ':' num2str(temp_second) ' min:sec' ])
else
temp_hour=fix(temp/3600);
temp=temp-temp_hour*3600;
temp_minute=fix(temp/60);
temp_second=mod(temp,60);
disp([' finish run with ' num2str(temp_hour) ':' num2str(temp_minute) ':' num2str(temp_second) ' h:m:sec' ])
end %temp<=60
end %for delay_i=1:delay_num
end %for fre=1000:1000:10000
case 'Save Threshold'
if 1
blockth.Filename=[zeng2.file(1:(length(zeng2.file)-2)) '_Z-' num2str(zeng.Iexstim(1).xyz(3)) '_Fre-' num2str(blockth.freq) '_Delay' num2str(zeng.Iinstim.delay) '_dt-' num2str(zeng.dt) '_D2-20.thr'];
blockth.Path=[zeng2.cd '\'];
else
[Filename Path]=uiputfile([thr.cd '\*.thr']);
end
if ~strcmp('0',num2str(blockth.Filename)) & ~strcmp('0',num2str(blockth.Path));
Have_dot=find(blockth.Filename=='.');
if isempty(Have_dot)
blockth.Filename=[blockth.Filename '.thr'];
end
dlmwrite([blockth.Path blockth.Filename],'');
FID=fopen([blockth.Path blockth.Filename],'w');
fprintf(FID,['\n X \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(1)) '\t' ]);
end
fprintf(FID,['\n Y \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(2)) '\t' ]);
end
fprintf(FID,['\n Z \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(3)) '\t' ]);
end
fprintf(FID,['\n Time step (us) \t']);
fprintf(FID,[num2str(zeng.dt)]);
fprintf(FID,['\n I %err \t']);
fprintf(FID,[num2str(blockth.cal.Ierrpercent)]);
fprintf(FID,['\n Dia \t']);
for i=1:blockth.cal.Dia_length
fprintf(FID,[num2str(blockth.cal.Dia(i)) '\t' ]);
end
fprintf(FID,['\t']);
fprintf(FID,['\t' 'X']);
for i=1:blockth.cal.Dia_length
fprintf(FID,[ '\t' ]);
end
fprintf(FID,['\t' 'Z']);
for freq_i=1:blockth.freq_num;
fprintf(FID,['\n \n Delay-Fre \t']);
fprintf(FID,[num2str(blockth.freq_all(freq_i)) '\t' ]);
for delay_i=1:blockth.delay_num
fprintf(FID,['\n ']);
fprintf(FID,[num2str(blockth.delay_all(freq_i,delay_i)) '\t' ]);
for j=1:blockth.cal.X_number
for i=1:blockth.cal.Dia_length
fprintf(FID,[num2str(blockth.cal.Ith(freq_i,delay_i,j,i)) '\t' ]);
end
fprintf(FID,['\t']);
for i=1:blockth.cal.Dia_length
fprintf(FID,[num2str(blockth.cal.X(j,i)) '\t' ]);
end
fprintf(FID,['\t']);
if isempty(blockth.cal.Z)
else
for i=1:blockth.cal.Dia_length
fprintf(FID,[num2str(blockth.cal.Z(j,i)) '\t' ]);
end
end
end
end
end
Status=fclose(FID);
end
end
