function HodgkinHuxley
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% Hodgkin-Huxley model neuron: Response of Squid Axon to Current Injection
% Please see class handouts for more details about the different
% equations used in this simulation code.
%
% G. Blohm (2010)
%-------------------------------------------------------------------------

clear

%% INITIALIZATION
DT = 0.025;     % time step (ms)
TMAX = 60;      % simulation time (ms)
t = 0:DT:TMAX;

VREST = 0;      % resting potential (Note: the voltage scale is shifted so 
                % that the resting potential vanishes)
                % for realistic simulations the voltage scale has to be
                % shifted by -65mV (i.e. shift ENA, EK, ELEAK and VREST)

GNA = 120.0;	% NA channel conductivity (mS/cm^2)
GK = 36.0;	    % K channel conductivity (mS/cm^2)
GLEAK = 0.3;	% leakage "channel" conductivity (mS/cm^2)

ENA = 115.0;	% Na reversal potential (mV)
EK = -12.0;   	% K reversal potential (mV)
ELEAK = 10.613;	% leakage current reversal potential (mV)

C_M = 1.0;	    % membrane capacitance (uF/cm^2)

v = zeros(size(t)); % memory pre-allocation (for faster simulations)
m = zeros(size(t));
n = zeros(size(t));
h = zeros(size(t));

v(1) = VREST;                                           % initial resting potential v
m(1) = alpha_m(VREST)/(alpha_m(VREST)+beta_m(VREST));	% initial gating variable m
h(1) = alpha_h(VREST)/(alpha_h(VREST)+beta_h(VREST));	% initial gating variable h
n(1) = alpha_n(VREST)/(alpha_n(VREST)+beta_n(VREST));	% initial hating variable n

%% stimulation stregth
stim = 10;              % strength of stimulation current (muA/cm^2)
I = stim*(t>5 & t<=30); % time dependent stimulation array (between 5 and 30 ms)

%% main simulation loop STARTS
% step-wise integration of the Hodgkin-Huxley equations
for i=2:length(t)
    gNa = GNA * m(i-1)^3 * h(i-1); % Na channel conductance (voltage and history/time dependent)
    gK  = GK  * n(i-1)^4;          % K channel conductance (voltage and history/time dependent)
    % note: the leakage conductance is constant

    mdot = alpha_m(v(i-1))*(1-m(i-1)) - beta_m(v(i-1))*m(i-1); % differential equations governing the gating variables
    hdot = alpha_h(v(i-1))*(1-h(i-1)) - beta_h(v(i-1))*h(i-1);
    ndot = alpha_n(v(i-1))*(1-n(i-1)) - beta_n(v(i-1))*n(i-1);
    vdot = (I(i-1) - gNa*(v(i-1)-ENA) - gK*(v(i-1)-EK) - GLEAK*(v(i-1)-ELEAK))/C_M; % differential equation governing the axon potential

    m(i) = m(i-1) + mdot*DT;		% Euler integration (see tutorial on numerical integration online)
    h(i) = h(i-1) + hdot*DT;
    n(i) = n(i-1) + ndot*DT;
    v(i) = v(i-1) + vdot*DT;
end
%% main simulation loop ENDS

%% plot results
plot(t,v,'b-'); hold on;
plot(t,I,'r:')
str = sprintf('Inject = %2.2f \\muA/cm^2',stim);
legend('axon potential','stimulation')
title(str);
xlabel('t (ms)'); ylabel('V (mV)')
axis([0 TMAX min(v) max(v)])


%% ========================================================================
%% EMPIRICAL SUBFUNCTIONS
% they contribute to the gaiting control of the different channels

%% alpha_h
function rate = alpha_h(v)
rate = 0.07*exp(-v/20.0);

%% alpha_m
function rate = alpha_m(v)
eps = 1e-12;

rate = zeros(size(v));		% initialization

idx = find(abs(v-25) > eps);		% normal case
if(~isempty(idx))
    rate(idx) = 0.1*(-v(idx)+25.0) ./ (exp((-v(idx)+25.0)/10.0)-1.0);
end

idx = find(abs(v-25) <= eps);		% handle sigularity V=25.0
if(~isempty(idx))
    rate(idx) = 1.0;
end

%% alpha_n
function rate = alpha_n(v)
eps = 1e-12;

rate = zeros(size(v));

idx = find(abs(v-10) > eps);
if(~isempty(idx))
    rate(idx) = 0.01*(-v(idx)+10.0) ./ (exp((-v(idx)+10.0)/10.0) - 1.0);
end

idx = find(abs(v-10) <= eps);		% HANDLE SIGULARITY AT V=10.0
if(~isempty(idx))
    rate(idx) = 0.1;
end

%% beta_h
function rate = beta_h(v)
rate =  1.0 ./ (exp((-v+30.0)/10.0) + 1.0);

%% beta_m
function rate =  beta_m(v)
rate = 4.0*exp(-v/18.0);

%% beta_n
function rate = beta_n(v)
rate = 0.125*exp(-v/80.0);
