% This is a simple ode integrator for the HH equations
% it is set up to be easily modified to use for other dynamics
function hh
 
global gnabar gkbar Vna Vk Vl Inp Veq Tfact Istim tstim BCL

tstim = 0.005;  %This is the length of time of an initial stimulus

Veq = 0; %is the original HH formulation
%Veq = -65; % is the physiological formulation

Tfact = 1; % this is the correction factor for temperatures other that 6.3C
%Tfact = 0.977 % correspnds to 0C

%Tfact = 1.085 % corresponds to 30 C.

gnabar = 120.;
gkbar = 36;
Vna = Tfact*(115+Veq) - Veq;
Vk = Tfact*(-12+Veq)-Veq;
Vl = Tfact*(10.5988+Veq)-Veq;

tstep = .1;
rn0 = 0.8;

Inp = 3.47;  % this is the amplitude of the initial input stimulus
%Inp = 3.5;
%Inp = 13.4;
%Inp=0;
Istim =0;  %this is the amplitude of the stedy current stimulus

BCL = 20;

tstep = .01;
t_end = 20;
V = 6.5098;  %The voltage threshold
V = Veq;
m = 0.0529;
h = 0.5961;
n = 0.3177;

%specify the output points
tspan = [0:tstep:t_end];

s0 = [V;m;h;n];

[T,S] = ode23(@deRHS,tspan, s0, odeset('maxstep',1));  
V = S(:,1);
m = S(:,2);
h = S(:,3);
n = S(:,4);

S(end,3)
S(end,4)
 
figure(4)
plot(n,h)
figure(1)
plot(T,V,'linewidth',2)
axis([0 20 -20 120])
xlabel('time (ms)','fontsize',20)
ylabel('Potential (mV)','fontsize',20)
figure(2)
plot(T,m,'b',T,n,'r',T,h,'g', 'linewidth',2)
xlabel('time (ms)','fontsize',20)
ylabel('Gating variables','fontsize',20)
axis([0 20 0 1])
text(4.5,0.85,'m(t)','fontsize',18)
text(14,0.36,'n(t)','fontsize',18)

text(14,0.64,'h(t)','fontsize',18)



end

%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gnabar gkbar Vna Vk Vl Inp Veq   Istim tstim BCL

%the original HH formulation has 
V = s(1)-Veq;
m = s(2);
h = s(3);
n = s(4);

 
Input = Inp/cosh(mod(t,BCL)-BCL/10) +Istim
 
 
AM=.1*(25.-V)./(exp(.1*(25.-V))-1.);
BM=4.*exp(-V/18);
AH=.07*exp(-V/20.);
BH=1./(exp(.1*(30.-V))+1.);
AN=.01*(10.-V)./(exp(.1*(10.-V))-1.);
BN=.125*exp(-V/80.);

Minf = AM./(AM+BM);

Fm=AM*(1.-m)-BM*m;
Fh=AH*(1.-h)-BH*h;
Fn=AN*(1.-n)-BN*n;
gna=gnabar*m.^3.*h;
%gnainf=gnabar*Minf.^3.*h;
gk=gkbar*n.^4;
Icl=.3*(V-Vl);

Fv = -gna*(V-Vna)-gk*(V-Vk)-Icl;

s_prime = [Fv+Input Fm Fh Fn]';

end
