%
%  Model B -- Shilnikov Attractor
%

clear all
global G1 G2 G_Cl1 G_Cl2 V1 V2 r1 r2 i1 i2 E_Na E_K E_Cl R el C Iin eps


G1=1; G2=-0.22; V1=-.8; V2=-.1; E_Na = 0.6; 

r1=6.6; r2=-.636; i1=.14; i2=.52; E_K = - 0.7; 

G_Cl1=.0000; G_Cl2=0; E_Cl = 0; 

el=.05; C=.5;R=.18; 

Iin=0; eps=.0005;


T=0; % transcient time
Tend=800;% Numerical solutions
%Tend=100;
options=odeset('e','off','BDF','on','AbsTol',1e-8,'RelTol',1e-8);

xinit=[-0.0211   -0.1082    2.5004    0.1081];
[t,x]=ode15s('SB2',[0,Tend],xinit, options);

N=find(t>T); % Skip the first T interval solution 

%xinit=x(end,:)

figure(3)
hold off
%nullclines; 
plot3(x(N:end,2),x(N:end,1),x(N:end,4),'r-')
axis([-.2 0 -.3 .3  0 .6])
grid on
xlabel('I_A')
ylabel('V_C')
zlabel('I_{K,p}')

return

%IV-Characteristics of Na and K tt
figure(1)
hold off
V=-1.2:.005:2;
i=-.5:.005:1.5;
plot(V,Na(V-E_Na),'r-',K(i)+E_K,i,'b-',V,Cl(V-E_Cl),'g-',...
    [-3 3], [0 0], [0 0], [-3 3])
axis([-2 2 -0.5 1])

%return


figure(2)
hold off
N=find(t>T); % Skip the first T interval solution 
plot(t(N:end),x(N:end,1),'r-',t(N:end),x(N:end,2),'b-',t(N:end),x(N:end,3),'m-',...
    t(N:end),(x(N:end,2)+x(N:end,3))/2,'k-.',t(N:end),(x(N:end,3)-x(N:end,2))/2,'k--')
hold on
legend('V_C','I_A','I_S','A_{Na}','A_K')

plot(t(N:end),x(N:end,2),'b-',t(N:end),x(N:end,3),'m-',...
    'linewidth',2)
hold on
xlabel('Time','fontsize',16)

% Below adds solution curve to IV-Characteristics tt
%figure(1)
%hold on
%plot(x(N:end,1),x(N:end,4),'g-')

% Nullcline (sum of two IV-Characteristics and constant Iin) on the I_A V_C
% axis tt
figure(102)
hold off
V=-1.2:.005:2;
i=-1.2:.005:2;
plot(V,-Na(V-E_Na)-(V-E_K)/r1+Iin,'g-',.1*V,.1*V/R,'m-',...
    [-3 3], [0 0], [0 0], [-3 3])
axis([-.5 .6 -0.5 .1])

% Below adds solution curve to IV-Characteristics tt
figure(102)
hold on
plot(x(N:end,1),x(N:end,2),'r-')

%figure(103)
%hold on
%plot(x(N:end,1),x(N:end,4),'g-')

%return

%3D picture tt

%figure(4)
%hold off

%nullclineszoomin; plot3(x(N:end,2),x(N:end,1),x(N:end,4),'r-')

return

%%%%%%%%%%%%%%%%%%%%%

figure(5)
hold off

xinit=[-0.1123    0.0041    3.9999    0.0979];

I_S0=.5;I_S1=4;I_Sinc=(I_S1-I_S0)/3;
for I_S=I_S0:I_Sinc:I_S1
xinit(3)=I_S; 
[t,x]=ode15s('SB2',[0,Tend],xinit);
%[t,x,a,b,c]=ode15s('SB2',[0,Tend],xinit,options);

N=find(t>T);
figure(5)
hold on
plot3(x(N:end,3),x(N:end,1),x(N:end,2),'r-',[xinit(3)], [xinit(1)],   [xinit(2)],'r.','markersize',20)
%xinit=x(end,:);

end

grid on


xlabel('I_S','fontsize',16)
ylabel('V_C','fontsize',16)
zlabel('I_A','fontsize',16)
axis([I_S0-.2 I_S1+.2 -.15 .44 -.01 0.025])

return



%%%%%%%% Burst frequency and Isospike keyboard plot %%%%%%%%%%%%%

figure(6)
hold off

xinit=[-0.1123    0.0041    3.9999    0.0979];

burstfreq=[];
bursttorefractory=[];
isospike=[];
I_S0=.75;I_S1=4;I_Sinc=(I_S1-I_S0)/40;

for I_S=I_S0:I_Sinc:I_S1
    xinit(3)=I_S; 
    [t,x]=ode15s('SB2',[0,Tend+60],xinit,options);
    xinit=x(end,:);
    xinit(1)=xinit(1)-.0002;

[t,x,a,b,c]=ode15s('SB2',[0,Tend+60],xinit,options);

    %freq=[freq, 1/(a(2)-a(1))];
    %spiketorefractory=[spiketorefractory, (a(2)-a(1))/a(1)];

    xinit=x(end,:);
    xinit(1)=xinit(1)-.0002;
    
    bstart=find(c==1);
    bend=find(c==3);
    
    burstfreq=[burstfreq,1/(a(bend)-a(bstart))];
    bursttorefractory=[bursttorefractory, (a(bend)-a(bstart))/a(bstart)];
    
    spks=(c==2);
    isospike=[isospike, sum(spks)];
end

I_S=I_S0:I_Sinc:I_S1;

% multiple axis plot

plot(I_S,burstfreq,'ro',I_S,bursttorefractory,'bd',I_S,isospike.*burstfreq,'m^','markersize',6)

%plot(I_S,isospike.*burstfreq,'m^','markersize',6)
grid on

legend('Burst Frequency','Burst:Refractory-Period','Location','N')

hold on
%plot(I_S,burstfreq,'r-',I_S,bursttorefractory,'b-')


xlabel('I_S','fontsize', 16)
ylabel('Burst-Spike Frequencies','fontsize', 16)

axis([0.5 1.1*I_S1 0 2.8])
box off

ax1 = gca;
set(ax1,'XColor','k','YColor','k')

ax2 = axes('Position',get(ax1,'Position'),...
           'XAxisLocation','top',...
           'YAxisLocation','right',...
           'Color','none',...
           'XColor','k','YColor','k');
hold on
maxspike=max(isospike);
plot(I_S,isospike,'ks','markersize',6,'Parent',ax2)

plot(I_S,isospike,'k-','Parent',ax2)
ylabel('Spike Per Burst','fontsize', 16)
axis([0.5 1.1*I_S1 0 8])

xlimits = get(ax1,'XLim');
ylimits = get(ax1,'YLim');
xinc = (xlimits(2)-xlimits(1))/8;
yinc = (ylimits(2)-ylimits(1))/8;
     
set(ax1,'YTick',[ylimits(1):yinc:ylimits(2)])        
       

return
