%
% Michael Schwemmer
% Math 227
% Hodgkin-Huxley Model
% ====================
% 01-29-2007
%
% simulates the response of the squid axon to an applied current.
%
% parameters as in 
%   Hoppensteadt and Peskin
%   Modeling and Simulation in Medicine and the Life Sciences
%
% to run, type the following at the matlab prompt: 
%               HH1(vstart,iappl0,istart0,istop0) 
%    vstart0 is the initial voltage, tmax0 is the total duration 
%    of the simulation, iappl0 is the amplitude of the applied 
%    current, istart is the time that the applied current is turned 
%    on, istop0 is the time that the applied current is turned off.
%
% voltages in mV, current in uA/cm^2, conductance in mS/cm^2, time is msec
%

function HH1(vstart0,tmax0,iappl0,istart0,istop0)

global C gNabar gKbar gLbar ENa EK EL istart istop iappl


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variables:
% membrance capacitance  uF/cm^2
C = 1.0; 
% Max Na conductance  mS/cm^2
gNabar = 120;
% Max K conductance  mS/cm^2
gKbar = 36;
% Max leakage conductance  mS/cm^2
gLbar = 0.3;
% Na Reversal Potential  mV
ENa = 45;
% K Reversal Potential  mV
EK = -82;
% Leakage Reversal Potential  mV
EL = -59;

vrest = -70;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% total simulations time  msec
tmax = tmax0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Current pulse parameters  uA/cm^2
iappl = iappl0; 
istart = istart0;
istop = istop0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial Values for m, n, and h are set to rest state

v = vrest;

m = alpham(v)/(alpham(v) + betam(v));
n = alphan(v)/(alphan(v) + betan(v));
h = alphah(v)/(alphah(v) + betah(v));

vstart=vstart0;
v=vstart;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the equations using ode45

% Initial Conditions 
s0 = [m n h v];

% Time span
tspan = [0,tmax];

options=odeset('InitialStep',10^(-3),'MaxStep',10^(-1));
[T,S] = ode45(@fxn,tspan,s0,options);

clf
figure(1)
subplot(2,1,1), plot(T,S(:,4)), xlabel('time'), ylabel('Voltage')
hold on
plot([istart,istop],[-90,-90],'r','linewidth',4)
subplot(2,1,2), plot(T,S(:,[1:3])), legend('m','n','h')



        function a=alpham(v)
            theta = (v+45)/10;
            if (theta ==0)
                a = 1.0;
            else
                a = 1.0*theta/(1-exp(-theta));
            end
        end
    
        function b = betam(v)
            b = 4.0*exp(-(v+70)/18);
        end
    
        function a = alphah(v)
            a = 0.07*exp(-(v+70)/20);
        end
    
        function b = betah(v)
            b = 1.0/(1+exp(-(v+40)/10));
        end
    
        function a=alphan(v)
            theta = (v+60)/10;
            if (theta ==0)
                a = 0.1;
            else
                a = 0.1*theta/(1-exp(-theta));
            end
        end

        function b = betan(v)
            b = 0.125*exp(-(v+70)/80);
        end
    
        
        function ds = fxn(t,s)
        
        global C gNabar gKbar gLbar ENa EK EL istart istop iappl

            ds = zeros(4,1);
            ds(1) = alpham(s(4))*(1-s(1))-betam(s(4))*s(1);
            ds(2) = alphan(s(4))*(1-s(2))-betan(s(4))*s(2);
            ds(3) = alphah(s(4))*(1-s(3))-betah(s(4))*s(3);
            gNa = gNabar*(s(1)^3)*s(3);
            gK = gKbar*(s(2)^4);
                if (istart<t)&&(t<istop)
                    ix = iappl;
                else
                    ix = 0;
                end
            ds(4) = -(1/C)*(gNa*(s(4)-ENa) + gK*(s(4)-EK) + gLbar*(s(4)-EL)) + ix/C;
           
        end
    
    
end