function hodgkin_huxley

% solves the Hodgkin-Huxley equations using forward Euler

global model

dt=.01;
tmax=200;
% need integer ratio tmax/dt
nsteps=tmax/dt;
if (floor(nsteps)~=nsteps)
  fprintf(' not an integer number of steps nsteps = %g \n',nsteps)
end

% model: 1=HH 2=connor-stevens (CoWa77)
model=1;

Iinj0=0;

time1=10;
duration1=200;
Iinj1=6.2;

time2=225;
duration2=5;
Iinj2=0;

Iinj=Iinj0;

% pre-allocate memory
 maxnumsteps=100000;
 t(1:maxnumsteps)=dt;
 Y(1:maxnumsteps,6)=0; 

% rest state
 V(1)=-65;
 n(1)=.32;
 m(1)=0.053;
 h(1)=.60;
 p(1)=0;
 q(1)=1;
 
 Y(1,:)=[V(1),n(1),m(1),h(1),p(1),q(1)];
 
 k=1;
 
 while t(k) < tmax
  
  if t(k) >= time1
    Iinj=Iinj1;
  end 
  if t(k) >= time1+duration1
    Iinj=Iinj0;
  end 
  if t(k) >= time2
    Iinj=Iinj2;
  end 
  if t(k) >= time2+duration2
    Iinj=Iinj0;
  end

  [Y(k+1,:)]=FE(Y(k,:),dt,Iinj);
  t(k+1)=t(k)+dt;
  k=k+1;

  if mod(k,1000) == 0
    fprintf('%g steps t = %g dt =%g \n',k,t(k),dt);
  end
  
 end   
 
numsteps=k-1;

% determine periods
V=Y(:,1);
Vsignchange=V(1:end-1).*V(2:end);
indsignchange=find(Vsignchange<0);
periods=indsignchange(3:end)-indsignchange(1:end-2);
periods=dt*periods

% plotting
  figure(1);
  subplot(3,1,1)
  plot(t(1:numsteps),Y(1:numsteps,1));
  axis([0 tmax -80 50]);
  legend(' voltage');
  subplot(3,1,2)
  plot(t(1:numsteps),Y(1:numsteps,2),'r',t(1:numsteps),Y(1:numsteps,3),'b',...
       t(1:numsteps),Y(1:numsteps,4),'m');
  legend('n','m','h');
  axis([0 tmax 0 1]);
  subplot(3,1,3)
  plot(t(1:numsteps),Y(1:numsteps,5),'k',t(1:numsteps),Y(1:numsteps,6),'g');  
  xlabel('Time')
  axis([0 tmax 0 1]);
  legend('p','q');
  
end

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

function Yp=F(Y,Iinj)

global model

if (model==1)
  Yp=F_hh(Y,Iinj);
else
  Yp=F_connor_stevens(Y,Iinj);
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Yp = F_hh(Y,Iinj) 
 
V=Y(1);
n=Y(2);
m=Y(3);
h=Y(4);
p=Y(5);
q=Y(6);

Gna=120;
Ena=50;

Gk=36;
Ek=-77;

EA=-75;
GA=0;

Cm=1;
Gl=.3;
El=-54.387;

Vthresh=-40;
alpham=0.1*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betam=4*exp(-0.0556*(V+65));
mp=alpham*(1-m)-betam*m;
 
alphah=0.07*exp(-0.05*(V+65));
betah=1/(1+exp(-0.1*(V+35)));
hp=alphah*(1-h)-betah*h;

Vthresh=-55;
alphan=0.01*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betan=0.125*exp(-0.0125*(V+65));
ninfinity=alphan/(alphan+betan);
taun=1/(alphan+betan);
np=(ninfinity-n)/taun;

alphap=2.5/(1+exp((V+20)*(-0.05)));
betap=0.5/(1+exp((V+0)*(0.08)));
pinfinity=alphap/(alphap+betap);
taup=1/(alphap+betap);
pinfinity=(0.0761*exp(0.0314*(V+94.22))/(1+exp(0.0346*(V+1.17))))^(1/3);
taup=0.3632+1.158/(1+exp(0.0497*(V+55.96)));
pp=(pinfinity-p)/taup;

alphaq=0.2/(1+exp((V+80)*(0.1)));
betaq=0.8/(1+exp((V+50)*(-0.05)));
qinfinity=alphaq/(alphaq+betaq);
tauq=1/(alphaq+betaq);
qinfinity=(1/(1+exp(0.0688*(V+53.3))))^4;
tauq=1.24+2.678/(1+exp(0.0624*(V+50)));
qp=(qinfinity-q)/tauq;

Vp=1/Cm*(Gna*m^3*h*(Ena-V)+Gk*n^4*(Ek-V)+Gl*(El-V)+Iinj+GA*p^3*q*(EA-V));

Yp(1)=Vp;
Yp(2)=np;
Yp(3)=mp;
Yp(4)=hp;
Yp(5)=pp;
Yp(6)=qp;

end

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

function Yp = F_connor_stevens(Y,Iinj) 
 
V=Y(1);
n=Y(2);
m=Y(3);
h=Y(4);
p=Y(5);
q=Y(6);

Gna=120;
Ena=55;

Gk=20;
Ek=-72.1;

EA=-75;
GA=47;

Cm=1;
Gl=.3;
El=-17;

Vthresh=-29.7;
alpham=0.38*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betam=15.2*exp(-0.0556*(V+54.7));
mp=alpham*(1-m)-betam*m;
 
alphah=0.266*exp(-0.05*(V+48));
betah=3.8/(1+exp(-0.1*(V+18)));
hp=alphah*(1-h)-betah*h;

Vthresh=-45.7;
alphan=0.02*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betan=0.25*exp(-0.0125*(V+55.7));
ninfinity=alphan/(alphan+betan);
taun=1/(alphan+betan);
np=(ninfinity-n)/taun;

alphap=2.5/(1+exp((V+20)*(-0.05)));
betap=0.5/(1+exp((V+0)*(0.08)));
pinfinity=alphap/(alphap+betap);
taup=1/(alphap+betap);
pinfinity=(0.0761*exp(0.0314*(V+94.22))/(1+exp(0.0346*(V+1.17))))^(1/3);
taup=0.3632+1.158/(1+exp(0.0497*(V+55.96)));
pp=(pinfinity-p)/taup;

alphaq=0.2/(1+exp((V+80)*(0.1)));
betaq=0.8/(1+exp((V+50)*(-0.05)));
qinfinity=alphaq/(alphaq+betaq);
tauq=1/(alphaq+betaq);
qinfinity=(1/(1+exp(0.0688*(V+53.3))))^4;
tauq=1.24+2.678/(1+exp(0.0624*(V+50)));
qp=(qinfinity-q)/tauq;

Vp=1/Cm*(Gna*m^3*h*(Ena-V)+Gk*n^4*(Ek-V)+Gl*(El-V)+Iinj+GA*p^3*q*(EA-V));

Yp(1)=Vp;
Yp(2)=np;
Yp(3)=mp;
Yp(4)=hp;
Yp(5)=pp;
Yp(6)=qp;

end

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

function [Y]=FE(Y,dt,Iinj)

global model

Y=Y+dt*F(Y,Iinj);

end

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