function varargout = FEA_Ve_HFBlock(varargin)
%find the high frequency blocking threshold for the below model file.
%--------------------------------------------------------------------------
%create('axon');
%global axon
%axon.nseg =151; %(15mm for 1um axon)
%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 HFblock
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=2; %1=calculated, 2=imported
zexst('change options','vestim',zeng2.options.vestim);
%------------------------------------
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=1; %us
zeng.tf=30; %ms
zeng2.rundisp =0;
%------------------------------------
Iexstim=zeng.Iexstim;
Iinstim=zeng.Iinstim;
%------------------------------------
FEAVe=get(findobj('type','figure','tag','ZengAnalysis','Name','FEA_Ve'),'userdata');
HFblock.feafile=FEAVe.Group.Filename;
%------------------------------------
HFblock.detect_start = 10;%ms
if HFblock.detect_start > Iinstim.delay
zexst('err',['The detect-start time takes place too early.'])
return
end
HFblock.detect_start_i = round(HFblock.detect_start/zeng.dt*1000);
HFblock.cal.Ierrpercent = 0.5; %the maximum percent err for the calculated threshold
HFblock.freq_all = [2000:2000:20000];
%HFblock.freq_all = [5000];
HFblock.freq_num = length(HFblock.freq_all);
HFblock.delay_num = 1;
HFblock.delay_all = zeros(HFblock.freq_num,HFblock.delay_num);
%------------------------------------
if 0
HFblock.cal.dia=10;
else
% HFblock.cal.dia=2:2:20;
HFblock.cal.dia=[5 15];
end
HFblock.cal.dia_length=length(HFblock.cal.dia);
%-----------------------------------------
%Istart at different value.
HFblock.Istart = 30000;
HFblock.Ilimit = 100000;%100mA
%-----------------------------------------
HFblock.cal.x_number=10;
HFblock.cal.x =zeros(HFblock.cal.x_number,HFblock.cal.dia_length);
HFblock.cal.yz=zeros(HFblock.cal.x_number,HFblock.cal.dia_length,2);
if 0 %systematic 0-L
HFblock.cal.x(:,1)=linspace(0,1,HFblock.cal.x_number)';
if HFblock.cal.dia_length>1
for i=2:HFblock.cal.dia_length
HFblock.cal.x(:,i)=HFblock.cal.x(:,1);
end
end
elseif 0 %systematic 0-L/2
HFblock.cal.x(:,1)=linspace(0,0.5,HFblock.cal.x_number)';
if HFblock.cal.dia_length>1
for i=2:HFblock.cal.dia_length
HFblock.cal.x(:,i)=HFblock.cal.x(:,1);
end
end
elseif 1 %random
HFblock.cal.x=rand(HFblock.cal.x_number,HFblock.cal.dia_length)*2-1; %[-1 +1]
end
eval(['Linter_D=' zeng2.var{1}.name '.Linter_D;'])
for i=1:HFblock.cal.dia_length
HFblock.cal.x(:,i)=HFblock.cal.x(:,i)*Linter_D*HFblock.cal.dia(i);
end
tic
if 0 %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
%-----------------------------------------
HFblock.cal.Ith=zeros(HFblock.freq_num,HFblock.delay_num,HFblock.cal.x_number,HFblock.cal.dia_length);
Total_axons=num2str(HFblock.cal.dia_length*HFblock.cal.x_number*HFblock.freq_num*HFblock.delay_num);
Total_axons_i=0;
%------------------------------------
HFblock.freq_i =1;
HFblock.cal.x_i =1;
HFblock.cal.dia_i =1;
disp(['---------------------------'])
disp(['dt = ' num2str(zeng.dt) 'us'])
disp(['tf = ' num2str(zeng.tf) 'ms'])
for x_i=1:HFblock.cal.x_number
HFblock.cal.x_i=x_i;
for freq_i=1:HFblock.freq_num;
HFblock.freq_i=freq_i;
HFblock.freq = HFblock.freq_all(freq_i); %this one is used in Isinth
HFblock.delay_all(freq_i,:)=Iinstim.delay+[0:(HFblock.delay_num-1)]/HFblock.delay_num*1000/HFblock.freq;
for dia_i=1:HFblock.cal.dia_length
HFblock.cal.dia_i=dia_i;
%---------------------------------------
%Assign a new axon position
FEAVe=get(findobj('type','figure','tag','ZengAnalysis','Name','FEA_Ve'),'userdata');
set(FEAVe.edit.x,'string',num2str(HFblock.cal.x(x_i,i)));
if 1
%random
disp(['Axon is randomly located within the fascicle.'])
FEAVe.Axon_i=[1 1 ceil(FEAVe.setting.number_Sites*rand)];
else
%at a fixed location
disp(['Axon is fixed at a location.'])
FEAVe.Axon_i=[1 1 find(FEAVe.setting.yz(:,1)==500 & FEAVe.setting.yz(:,2)==0)];
end
set(FEAVe.figure,'userdata',FEAVe)
HFblock.cal.yz(x_i,dia_i,:)=FEAVe.setting.yz(FEAVe.Axon_i(3),:);
%-----------------------------------------
eval([zeng2.var{1}.name '.dia =HFblock.cal.dia(dia_i);'])
if 1 %Fix L
temp.L=51000; %um
eval([zeng2.var{1}.name '.nseg =1+2*round(temp.L/2/HFblock.cal.dia(dia_i)/' zeng2.var{1}.name '.Linter_D);']);
eval([zeng2.var{1}.name '.L =' zeng2.var{1}.name '.nseg*HFblock.cal.dia(dia_i)*' zeng2.var{1}.name '.Linter_D;']);
disp(['L is fixed at ' num2str(temp.L) ' um'])
end
for delay_i=1:HFblock.delay_num
zeng.Iinstim.delay=HFblock.delay_all(freq_i,delay_i);
%-------------------------------------------------------------------------------------------------------------------------------
Total_axons_i=Total_axons_i+1;
disp([num2str(Total_axons_i) '/' Total_axons])
disp([' Fre=' num2str(HFblock.freq) ' delay=' num2str(zeng.Iinstim.delay) ' Dia=' num2str(HFblock.cal.dia(dia_i)) ' nseg=' num2str(axon.nseg) ' x=' num2str(HFblock.cal.x(x_i,dia_i)) ' yz=' num2str(FEAVe.setting.yz(FEAVe.Axon_i(3),:))])
if Total_axons_i==1
Iminmaxmean=HFblock.Istart;
else
%------For freq_i=1----------------
%Iminmaxmean = the Ith of previous simulation
%It was automatically assigned at the end of the previous
%simulation
%-----For freq_i>1-----------------------------
%Iminmaxmean = Ith of previous frequency
if freq_i>1
Iminmaxmean=HFblock.cal.Ith(freq_i-1,delay_i,x_i,dia_i);
end
%-----------------------------------------
if Iminmaxmean>=HFblock.Ilimit
Iminmaxmean=HFblock.Ilimit/8;
end
end
Iminmax(1:2)=Iminmaxmean;
zeng2.amp =Iminmaxmean;
zexst('calculate');
if zeng2.notreadytorun
return
end
if zeng2.singular
Blocked=0;
Iminmax(2)=HFblock.Ilimit;
else
if isempty(find(axon.vm(1,[HFblock.detect_start_i:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
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)*1.25;
zeng2.amp =Iminmax(2);
zexst('calculate');
if zeng2.singular
Blocked=0;
Iminmax(2)=HFblock.Ilimit;
else
if isempty(find(axon.vm(1,[HFblock.detect_start_i:size(axon.vm,2)])>0))
Blocked=1;
else
Blocked=0;
end
end
disp([' I = ' num2str([zeng2.amp]) ' - ' num2str(Blocked)])
if Iminmax(2)>=HFblock.Ilimit
disp([' I > the limit (' num2str(HFblock.Ilimit) ' uA)'])
%These two lines help skip the while loop below
Iminmax(1)=HFblock.Ilimit;
Iminmax(2)=HFblock.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,[HFblock.detect_start_i: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)>=HFblock.Ilimit
HFblock.cal.Ith(freq_i,delay_i,x_i,dia_i)=HFblock.Ilimit;
break
end
zeng2.amp =Iminmaxmean;
zexst('calculate');
if isempty(find(axon.vm(1,[HFblock.detect_start_i: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)>HFblock.cal.Ierrpercent
Iminmax(2)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmaxmean])])
HFblock.cal.Ith(freq_i,delay_i,x_i,dia_i)=Iminmaxmean;
break
end
elseif Blocked==0
if (100*(Iminmax(2)-Iminmaxmean)/Iminmaxmean)>HFblock.cal.Ierrpercent
Iminmax(1)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmax(2)])])
HFblock.cal.Ith(freq_i,delay_i,x_i,dia_i)=Iminmax(2);
break
end
end
end
FEA_Ve_HFBlock('Save Threshold')
end %for x_i=1:HFblock.cal.x_number
FEA_Ve_HFBlock('Save Threshold')
end% for i=1:HFblock.cal.dia_length
FEA_Ve_HFBlock('Save Threshold')
temp=toc;
if 0
figure
if size(HFblock.cal.Ith(freq_i,delay_i,:,:),1)==1
plot(HFblock.cal.dia,HFblock.cal.Ith(freq_i,delay_i,:,:));
else
errorbar(HFblock.cal.dia,mean(HFblock.cal.Ith(freq_i,delay_i,:,:)),std(HFblock.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
HFblock.Filename=[zeng2.file(1:(length(zeng2.file)-2)) '_sim-' num2str(HFblock.cal.x_i) '_fre-' num2str(HFblock.freq_all(HFblock.freq_i)) '_dia-' num2str(HFblock.cal.dia(HFblock.cal.dia_i)) '_dt-' num2str(zeng.dt) '_D1-20.thr'];
HFblock.Path=[zeng2.cd '\'];
else
[Filename Path]=uiputfile([thr.cd '\*.thr']);
end
if ~strcmp('0',num2str(HFblock.Filename)) & ~strcmp('0',num2str(HFblock.Path));
Have_dot=find(HFblock.Filename=='.');
if isempty(Have_dot)
HFblock.Filename=[HFblock.Filename '.thr'];
end
dlmwrite([HFblock.Path HFblock.Filename],'');
FID=fopen([HFblock.Path HFblock.Filename],'w');
fprintf(FID,['\n' HFblock.feafile]);
fprintf(FID,['\n' zeng2.file]);
fprintf(FID,['\n Time step (us)']);
fprintf(FID,['\t' num2str(zeng.dt)]);
fprintf(FID,['\n Ierr percent']);
fprintf(FID,['\t' num2str(HFblock.cal.Ierrpercent)]);
for x_i=1:HFblock.cal.x_i
fprintf(FID,['\n' '\n' num2str(x_i)]);
fprintf(FID,['\n Dia']);
for i=1:HFblock.cal.dia_length
fprintf(FID,['\t' num2str(HFblock.cal.dia(i)) ]);
end
fprintf(FID,['\t' 'X']);
for dia_i=1:HFblock.cal.dia_length
fprintf(FID,['\t' num2str(HFblock.cal.x(x_i,dia_i))]);
end
fprintf(FID,['\t' 'Y']);
for dia_i=1:HFblock.cal.dia_length
fprintf(FID,['\t' num2str(HFblock.cal.yz(x_i,dia_i,1))]);
end
fprintf(FID,['\t' 'Z']);
for dia_i=1:HFblock.cal.dia_length
fprintf(FID,['\t' num2str(HFblock.cal.yz(x_i,dia_i,2))]);
end
for freq_i=1:HFblock.freq_num;
fprintf(FID,['\n Delay-Fre']);
fprintf(FID,['\t' num2str(HFblock.freq_all(freq_i))]);
for delay_i=1:HFblock.delay_num
fprintf(FID,['\n' num2str(HFblock.delay_all(freq_i,delay_i))]);
for dia_i=1:HFblock.cal.dia_length
fprintf(FID,[ '\t' num2str(HFblock.cal.Ith(freq_i,delay_i,x_i,dia_i))]);
end
end
end
end
Status=fclose(FID);
end
end
