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

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

istart = 150; % time applied current starts
ilength=300;   % length of applied current pulse
Ie=0;     % magnitude of applied current pulse
I0= 10;     % magnitude of applied current outside pulse --- set to 10 to see "anode break"
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

t=0:dt:tmax;    % time vector
V=zeros(size(t)); % 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(size(t));   % n: potassium activation gating variable
n(1) = 0.0;         % start off at zero
m=zeros(size(t));   % m: sodium activation gating variable
m(1) = 0.0;         % start off at zero
h=zeros(size(t));   % h: sodim inactivation gating variable
h(1) = 0.0;         % start off at zero

Iapp=zeros(size(t)); % Applied current, relevant in current-clamp mode
I_L=zeros(size(t)); % in case we want to plot and look at the total current
I_Na=zeros(size(t)); % in case we want to plot and look at the total current
I_K=zeros(size(t)); % in case we want to plot and look at the total current

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

Vapp=zeros(size(t)); % Applied voltage, relevant in voltage-clamp mode
if ( vclamp_flag )
    for i = 1:round(vstart/dt)          % % make V0 before pulse
        Vapp(i) = V0;
    end
    for i=round(vstart/dt)+1:round((vstart+vlength)/dt) % make Ve for duration of voltage pulse
        Vapp(i) = Ve;
    end
    for i=round((vstart+vlength)/dt):length(Vapp) % make V0 following pulse
        Vapp(i) = V0;
    end
end

Itot=zeros(size(t)); % in case we want to plot and look at the total current

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

    Vm = V(i-1);

    % 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(1) = m_inf;
        h(1) = h_inf;
        n(1) = n_inf;
    end
    
    m(i) = m(i-1) + (m_inf-m(i-1))*dt/tau_m;    % Update m

    h(i) = h(i-1) + (h_inf-h(i-1))*dt/tau_h;    % Update h

    n(i) = n(i-1) + (n_inf-n(i-1))*dt/tau_n;    % Update n

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

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

    Itot(i-1) = I_L(i-1)+I_Na(i-1)+I_K(i-1)+Iapp(i-1); % total current is sum of leak + active channels + applied current

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


    if ( vclamp_flag )      % if we are using voltage clamp
        V(i) = Vapp(i);     % ignore the voltage integration and set V to be the applied voltage
    end

    if ( spikenow == 0 )
        if ( V(i) < -60 )
            if ( ( t(i) > istart ) && ( t(i) < istart+ilength ) )
                spikes(i) = 1;
            end
            spikenow = 1;
        end
    else
        if ( V(i) > -50 )
            spikenow = 0;
        end
    end


end

hold on;

plot(t,V);