% Hodgkin_Huxley_model.m
% APMA 0410
% Fall 2008

% Named after Alan Hodgkin and Sir Andrew Huxley
%   (Sir Andrew is alive and well at 90,
%     he is the brother of Aldous Huxley, of Brave New World fame).

% H and H won the Nobel prize in medicine in 1963
% for their study (in the 30s and 40s) of the action potential
%    (in the frog's sciatic nerve and in the giant squid axon). 

% They also proposed the HH model,
%    a mathematical/computational model,
%    which LATER proved to be correct to a remarkable extent,
%    and which spawned a lot of theoretical work.

% http://en.wikipedia.org/wiki/Andrew_Huxley

% The HH model is a coupled system of 4 nonlinear differential equations
%    for the following four variables:

%  V(t) (membrane potential),
%  m(t) (gating variable),
%  h(t) (gating variable),
%  n(t) (gating variable).


% The differential equation for V(t) is the simplest of all possible ODE's:
%     1-st order and linear.
%     It models a single RC circuit.
%     (The HH model is therefore called a single-compartment model.)

% It has CONSTANT coefficients

%     EXCEPT FOR THE FACT THAT

% The current i(t) that is injected into the membrane is the sum of three
%   currents:

%   a leakage current (simple RC- behavior, i.e. constant coefficients)

%   a K+ current

%   an Na+ current

% The latter two depend on V itself in a highly NONLINEAR fashion,
%    through the gating variables m, h, and n.

%  HOW SO?

% Each one of the gating variables m, h and n obeys
%    a 1-st order linear equation,
%    whose coefficients are non-linear functions of V itself.

%  H and H PREDICTED the shape and parameters of these non-linear functions
%    in a very accurate fashion.


% The code below numerically integrates this ENTIRE system
%    of coupled differential equations


% The code allows for three different types of experiment,
%    corresponding to different setups for the external current,
%         a.k.a. "electrode current":

%    EXPERIMENT 1: use a constant positive current to excite
%         (i.e., "depolarize") the neuron.

%    EXPERIMENT 2: use a brief pulse of negative current.

%    EXPERIMENT 3: use a random time-varying current.


clear
close all

%  PARAMETERS OF MODEL

epoch = .50;                        % length of epoch in seconds
bin_length = .0001;                 % bin length (discretization of time axis) in sec

E_L  = -.054387;                    % leakage reversal potential (V)
E_K  = -.077;                       % K+ reversal potential (V)
E_Na =  .050;                       % Na+ reversal potential (V)
g_bar_L  = 3*10^-6;                 % leakage conductance (S/mm^2)
g_bar_K  = .36*10^-3;               % maximal K+ conductance (S/mm^2)
g_bar_Na = 1.2*10^-3;               % maximal Na+ conductance (S/mm^2)
c_m = 10^-8;                        % specific membrane capacitance (F/mm^2)
                                    %       (note: units are consistent, because of the use of mm^2 for area)

time_axis = 0:bin_length:epoch;     % time axis
N = length(time_axis);              % number of bins in spike train

current_setup = 1;                  % see below

switch current_setup
    
case 1                                  % constant electrode current, ranging from 0 to 500 nA/mm^2
    
    Ie_values = [0 50 60 62 63 65 100:50:500]*10^-9;       % electrode current values (A/mm^2)
    window = .005/bin_length:N;         % current-injection window in bin units
    window = round(window);
    
case 2                                  % pulse of negative current followed by zero current
    
    Ie_values = -50*10^-9;              % electrode current value  (A/mm^2)
    window = 1:.025/bin_length;         % current-injection window in bin units
    
case 3                                  % white noise current
    
    Ie_values = 1;                      % dummy
    sampling_f = 1/bin_length;          % sampling frequency in Hz
    delta_t = .0003;                    % time-step interval in sec
    sigma_s = 4*10^-9;                  % standard deviation of current in A/mm^2
    m = delta_t*sampling_f;             % number of sampling points in each time-step interval
    m = round(m);
    current = randn(1,round(N/m));                % i.i.d. standard gaussian rv's
    current = kron(current,ones(1,m));            % replicate each value m times
    current = current*sigma_s;                    % multiply by sigma_s
    current = current/sqrt(delta_t);              % divide by sqrt(delta_t) (see book p. 23)
    N = length(current);
    window = 1:N;                                 % current-injection window in bin units
    time_axis = window*sampling_f;
    
end

V = -.065*ones(1,N);
m = 0.0529*ones(1,N);
h = 0.5961*ones(1,N);
n = 0.3177*ones(1,N); 
i_m = zeros(1,N);

k = 0;
figure_number = 0;

for Ie = Ie_values
    
    k = k+1;
    
    if current_setup < 3        
        current = zeros(1,N);               % electrode current will be zero except in window
        current(window) = Ie;        
    end
    
    for b = 1:N-1
        
        v = V(b)*1000;                                          % V expressed in mV
        alpha_n = .01*(v+55)/(1-exp(-.1*(v+55)));               % eq. 5.22
        beta_n  = .125*exp(-.0125*(v+65));                      % eq. 5.22
        
        alpha_m =  .1*(v+40)/(1-exp(-.1*(v+40)));               % eq. 5.24
        beta_m  =   4*exp(-.0556*(v+65));                       % eq. 5.24
        
        alpha_h = .07*exp(-.05*(v+65));                         % eq. 5.24
        beta_h  = 1/(1+exp(-.1*(v+35)));                        % eq. 5.24
        
        tau_n = .001/(alpha_n+beta_n);                          % eq. 5.18 expressed in sec
        tau_m = .001/(alpha_m+beta_m);                          % eq. 5.18 expressed in sec
        tau_h = .001/(alpha_h+beta_h);                          % eq. 5.18 expressed in sec
        
        n_infinity = alpha_n/(alpha_n+beta_n);                  % eq. 5.19
        m_infinity = alpha_m/(alpha_m+beta_m);                  % eq. 5.19
        h_infinity = alpha_h/(alpha_h+beta_h);                  % eq. 5.19
        
        n(b+1) = n_infinity + (n(b) - n_infinity)*exp(-bin_length/tau_n);     % eq. 5.48
        m(b+1) = m_infinity + (m(b) - m_infinity)*exp(-bin_length/tau_m);     % eq. 5.48
        h(b+1) = h_infinity + (h(b) - h_infinity)*exp(-bin_length/tau_h);     % eq. 5.48
        
        i_m(b) = g_bar_L*(V(b)-E_L) + g_bar_K*(n(b)^4)*(V(b)-E_K) + g_bar_Na*(m(b)^3)*h(b)*(V(b)-E_Na);   % eq. 5.25
        
        sum_g = g_bar_L + g_bar_K*(n(b+1)^4) + g_bar_Na*(m(b+1)^3)*h(b+1);                     % sum of specific conductances 
        
        V_infinity = (g_bar_L*E_L + g_bar_K*(n(b+1)^4)*E_K + g_bar_Na*(m(b+1)^3)*h(b+1)*E_Na + current(b+1))/sum_g;   % eq. 5.49
        tau_V = c_m/sum_g;                                                                                            % eq. 5.50
        
        V(b+1) = V_infinity + (V(b) - V_infinity)*exp(-bin_length/tau_V);       % eq. 5.48
        
    end
    
    % spikes times are found by looking for a maximum of V at a positive value:
    spike_train = zeros(1,N);
    spike_train(1:N-2) = V(1:N-2)>0 & diff(V(1:N-1))>0 & diff(V(2:N))<0 ;
    % compute rate over current injection window:
    rate(k) = sum(spike_train(window))/(length(window)*bin_length); 
  
    if k < 6 || k ==15
        
        figure_number = figure_number + 1;
        figure(figure_number)
        subplot(5,1,1)
        [haxes,hline1,hline2] = plotyy(time_axis,V,time_axis,current);
        hold on, plot(time_axis,E_L*ones(1,N),'--',time_axis,E_K*ones(1,N),'--',time_axis,E_Na*ones(1,N),'--'), hold off    % plot the three reversal potentials
        axes(haxes(1))
        ylabel('membrane potential (V)')
        axes(haxes(2))
        ylabel('injected current (A/mm^2)')
        set(gca,'YLim',[1.2*min(current) max(eps,1.2*max(current))])
        subplot(5,1,2), plot(time_axis,i_m)
        grid
        ylabel('i_m')
        %set(gca,'YLim',[0 2.5]*10^-9)
        subplot(5,1,3), plot(time_axis,m)
        grid
        ylabel('m')
        %set(gca,'YLim',[0 2.5]*10^-9)
        subplot(5,1,4), plot(time_axis,h)
        grid
        ylabel('h')
        %set(gca,'YLim',[0 2.5]*10^-9)
        subplot(5,1,5), plot(time_axis,n)
        grid
        ylabel('n')
        %set(gca,'YLim',[0 2.5]*10^-9)
        xlabel('time (sec)')
        
    end
    
end

if current_setup ==1
    
    figure_number = figure_number + 1;
    figure(figure_number)
    plot(Ie_values,rate,'-*')
    title(['FIRING RATE v. INPUT CURRENT'])
    ylabel(['firing rate (Hz)'])
    xlabel('input current (A)')
    
end