%----------------------------------------------------------------------
%Adapted by Bruce Land from:
% Physiology 317 - Methods in Computational Neurobiology
% HODGKIN-HUXLEY: Response of Squid Axon to Current Injection
% M. Nelson  1995, 1997
%----------------------------------------------------------------------

clear
clf
orient landscape

%%%%%%%%%%%%%%%%%%%%%%%
% INITIALIZATION
%%%%%%%%%%%%%%%%%%%%%%%
DT = 0.025;
TMAX = 200;
t = 0:DT:TMAX;

VREST = -60.0;

GNA = 120.0;	% mS/cm^2
GK = 36.0;	% mS/cm^2
GLEAK = 0.3;	% mS/cm^2

ENA = 115.0 + VREST;	% mV
EK = -12.0 + VREST;	% mV
ELEAK = 10.613 + VREST;	% mV

C_M = 1.0;	% uF/cm^2

% %current strength
I = 10 ;
I2 = 10.1 ;

%synapse
GSYN = 5 ;
ESYN = EK ;

% %current strength 2:1
% I = 6.0 ;
% I2 = 2*I ;
% %synapse
% GSYN = 10 ;
% ESYN = EK ;

%%%%%%%%%%%%%%%%%%%%%%%
% MAIN SIMULATION LOOP
%%%%%%%%%%%%%%%%%%%%%%%

clear v  m  h  n				% clear old varibles
clear v2  m2  h2  n2				% clear old varibles
v = VREST;					% initial membrane voltage
v2 = VREST;					% initial membrane voltage

m = alpha_m(v-VREST)/(alpha_m(v-VREST)+beta_m(v-VREST));	% initial (steady-state) m
h = alpha_h(v-VREST)/(alpha_h(v-VREST)+beta_h(v-VREST));	% initial (steady-state) h
n = alpha_n(v-VREST)/(alpha_n(v-VREST)+beta_n(v-VREST));	% initial (steady-state) n

m2 = alpha_m(v-VREST)/(alpha_m(v-VREST)+beta_m(v-VREST));	% initial (steady-state) m
h2 = alpha_h(v-VREST)/(alpha_h(v-VREST)+beta_h(v-VREST));	% initial (steady-state) h
n2 = alpha_n(v-VREST)/(alpha_n(v-VREST)+beta_n(v-VREST));	% initial (steady-state) n

for i=2:length(t)

    M = m(i-1);	% get values from last time step
    H = h(i-1);	% (hopefully this substitution makes
    N = n(i-1);	% the following code a bit easier to read)
    V = v(i-1);

    M2 = m2(i-1);	% get values from last time step
    H2 = h2(i-1);	% (hopefully this substitution makes
    N2 = n2(i-1);	% the following code a bit easier to read)
    V2 = v2(i-1);

    gNa = GNA * M^3 * H;
    gK  = GK  * N^4;

    gNa2 = GNA * M2^3 * H2;
    gK2  = GK  * N2^4;

    gS  = GSYN * (V2>20) ;
    gS2 = GSYN * (V >20) ;

    mdot = alpha_m(V-VREST)*(1-M) - beta_m(V-VREST)*M;
    hdot = alpha_h(V-VREST)*(1-H) - beta_h(V-VREST)*H;
    ndot = alpha_n(V-VREST)*(1-N) - beta_n(V-VREST)*N;
    vdot = (I - gS*(V-ESYN)- gNa*(V-ENA) - gK*(V-EK) - GLEAK*(V-ELEAK))/C_M;

    mdot2 = alpha_m(V2-VREST)*(1-M2) - beta_m(V2-VREST)*M2;
    hdot2 = alpha_h(V2-VREST)*(1-H2) - beta_h(V2-VREST)*H2;
    ndot2 = alpha_n(V2-VREST)*(1-N2) - beta_n(V2-VREST)*N2;
    vdot2 = (I2 - gS2*(V2-ESYN) - gNa2*(V2-ENA) - gK2*(V2-EK) - GLEAK*(V2-ELEAK))/C_M;

    m(i) = m(i-1) + mdot*DT;		% Euler integration
    h(i) = h(i-1) + hdot*DT;
    n(i) = n(i-1) + ndot*DT;
    v(i) = v(i-1) + vdot*DT;

    m2(i) = m2(i-1) + mdot2*DT;		% Euler integration
    h2(i) = h2(i-1) + hdot2*DT;
    n2(i) = n2(i-1) + ndot2*DT;
    v2(i) = v2(i-1) + vdot2*DT;
end

figure(1)
plot(t,v,'b-'); hold on;
plot(t,v2,'r-'); hold on;

str = sprintf('Inject = %2.2f \\muA/cm^2',I);
title(str);
drawnow;