% This model is the Hodgkin-Huxley model of Homework 2.
clear
dt = 0.002;
tmax=60;

iclamp_flag = 1; % if this is 1, run under current clamp conditions
vclamp_flag = 0; % otherwise this should be 1, for voltage clamp conditions

xbins = 150;
istart = 10; % time applied current starts
ilength=40;   % length of applied current pulse
Ie= -200;     % magnitude of applied current pulse
I0= -20;     % magnitude of applied current outside pulse

V_L = -10.613;   % leak reversal potential
E_Na = -115;   % reversal for sodium channels
E_K = 12;   % reversal for potassium channels

g_L = 0.3;     % specific leak conductance
g_Na = 120;  % specific sodium conductance
g_K = 36;     % specific potassium conductance

cm = 1;     % specific membrane capacitance
tau = cm/g_L;   % membrane time constant

diffact = 100;    % factor proportional to conductance between neighboring sections of axon

t=0:dt:tmax;    % time vector
V=zeros(length(t),xbins); % voltage vector
spikes = zeros(size(t)); % will contain set of spike times
spikenow = 0;       % variable to say if code is in middle of a spike

if ( iclamp_flag ) % i.e. if in current-clamp mode
    V(1) = V_L;    % set the inititial value of voltage
end

n=zeros(1,xbins);   % n: potassium activation gating variable
n(1) = 0.0;         % start off at zero
m=zeros(1,xbins);   % m: sodium activation gating variable
m(1) = 0.0;         % start off at zero
h=zeros(1,xbins);   % h: sodim inactivation gating variable
h(1) = 0.0;         % start off at zero

Iapp=zeros(size(t)); % Applied current, relevant in current-clamp mode

if ( iclamp_flag )   % i.e. if in current-clamp mode
    for i=1:round((istart)/dt) % make non-zero for duration of current pulse
        Iapp(i) = I0;
    end
    for i=round(istart/dt)+1:round((istart+ilength)/dt) % make non-zero for duration of current pulse
        Iapp(i) = Ie;
    end
    for i=round((istart+ilength)/dt+1:length(Iapp)) % make non-zero for duration of current pulse
        Iapp(i) = I0;
    end
end

Itot=zeros(1,xbins); % in case we want to plot and look at the total current

for i = 2:length(t); % now see how things change through time
    for j = 1:xbins
        I_L = g_L*(V_L-V(i-1,j));

        Vm = V(i-1,j);

        if j == 1 
            d2vdx2 = V(i-1,2)-V(i-1,1);
        else
            if j == xbins
                d2vdx2 = V(i-1,xbins-1) - V(i-1,xbins);
            else
                d2vdx2 = V(i-1,j+1)+V(i-1,j-1)-2*V(i-1,j);
            end
        end
        
        Idiff = diffact*d2vdx2;
        
        % Sodium and potassium gating variables are defined by the
        % voltage-dependent transition rates between states, labeled alpha and
        % beta. Written out from Dayan/Abbott, units are 1/ms.

        if ( Vm == -25 )
            alpha_m = 0.1/0.1;
        else
            alpha_m = (0.1*(Vm+25))/(exp(0.1*(Vm+25))-1);
        end
        beta_m = 4*exp(Vm/18);
        alpha_h = 0.07*exp(Vm/20);
        beta_h = 1/(1+exp((Vm+30)/10));
        if ( Vm == -10)
            alpha_n = 0.01/0.1;
        else
            alpha_n = (0.01*(Vm+10))/(exp(0.1*(Vm+10))-1);
        end
        beta_n = 0.125*exp((Vm)/80);

        % From the alpha and beta for each gating variable we find the steady
        % state values (_inf) and the time constants (tau_) for each m,h and n.

        tau_m = 1/(alpha_m+beta_m);      % time constant converted from ms to sec
        m_inf = alpha_m/(alpha_m+beta_m);

        tau_h = 1/(alpha_h+beta_h);      % time constant converted from ms to sec
        h_inf = alpha_h/(alpha_h+beta_h);

        tau_n = 1/(alpha_n+beta_n);      % time constant converted from ms to sec
        n_inf = alpha_n/(alpha_n+beta_n);

        if ( i == 2 )
            m(j) = m_inf;
            h(j) = h_inf;
            n(j) = n_inf;
        end

        m(j) = m_inf - (m_inf-m(j))*exp(-dt/tau_m);    % Update m

        h(j) = h_inf - (h_inf-h(j))*exp(-dt/tau_h);    % Update h

        n(j) = n_inf - (n_inf-n(j))*exp(-dt/tau_n);    % Update n

        I_Na = g_Na*m(j)*m(j)*m(j)*h(j)*(E_Na-V(i-1,j)); % total sodium current

        I_K = g_K*n(j)*n(j)*n(j)*n(j)*(E_K-V(i-1,j)); % total potassium current

        % total current is sum of leak + active channels + applied current
        % + current diffusing from other positions on the axon
       Itot(j) = I_L+I_Na+I_K+Idiff; 
        if ( j == 1 )
            Itot(j) = Itot(j) + Iapp(i);
        end

        V(i,j) = V(i-1,j) + Itot(j)*dt/cm;         % Update the membrane potential, V.


    end
    if ( mod(i,100) == 1 )
        figure(1);
        axis([0 xbins -110 10]);
        plot(V(i,:))
        axis([0 xbins -110 10]);
        t(i)
    end
    
end

hold on;