Simulation of Hogkin-Huxley Neuron, Anode Break Potential

hh.m

close all;
clear all;
% feed positive current and check plots of v,m,h,n
[t , r] = ode15s(@hh,[0 10],[-60  0.052932 0.596121 0.317677]);

    figure
            subplot(2,1,1)
            plot(t,r(:,1));
            title('hogkin huxley')
            
            
            subplot(2,1,2)
             hold on
                plot(t,r(:,2));
                plot(t,r(:,3));
                plot(t,r(:,4));
                legend('m','h','n')
             hold off
    grid

% anode break potential - v,m,h,n
[t, r] = ode15s(@hh_i_negative_for_anode_break, [0 50], [-60  0.052932 0.596121 0.317677]);

    figure
            subplot(2,1,1)
            hold on
                plot(t,r(:,1));
                title('hogkin huxley - anode break')
                plot([0 10 10 20 30 30 40 50], [0 0 -20 -20 -20 0 0 0]);
            hold off
            
            subplot(2,1,2)
             hold on
                plot(t,r(:,2));
                plot(t,r(:,3));
                plot(t,r(:,4));
                legend('m','h','n')
             hold off
    grid


function result = hh(t,r)

    % vars
    g_k_bar = 36;
    e_k = -72;

    g_na_bar = 120;
    e_na = 55;

    g_l = 0.3;
    e_l = -49.401079;

    c = 1;
    iext = 10;

   

    if r(1) == -35
        alpha_m = 1;
    else
        alpha_m = (-0.1 * (r(1) + 35))/(exp(-(r(1) + 35)/10) - 1);
    end
    beta_m = 4 * exp(-(r(1) + 60)/18);


    if r(1) == -50 
        alpha_n = 0.1;
    else
        alpha_n = (-0.01 * (r(1) + 50))/(exp(-(r(1) + 50)/10) - 1);
    end
    beta_n = 0.125 * exp(-(r(1) + 60)/80);

    alpha_h = 0.07 * exp(-(r(1) + 60)/20);
    beta_h = 1/(1 + exp(-(r(1)+30)/10));
    
    result = zeros(4,1); % v,m,h,n
    result(1) = (1/c) * ( iext - (g_k_bar * r(4)^4 * (r(1) - e_k)) - (g_na_bar * r(2)^3 * r(3) * (r(1) - e_na)) - (g_l * (r(1) - e_l)) );
    result(2) = (alpha_m * (1 - r(2))) - (beta_m * r(2));   
    result(3) = (alpha_h * (1 - r(3))) - (beta_h * r(3));
    result(4) = (alpha_n * (1 - r(4))) - (beta_n * r(4));
end

function result = hh_i_negative_for_anode_break(t,r)

    % vars
    g_k_bar = 36;
    e_k = -72;

    g_na_bar = 120;
    e_na = 55;

    g_l = 0.3;
    e_l = -49.401079;

    c = 1;
   
   

    if r(1) == -35
        alpha_m = 1;
    else
        alpha_m = (-0.1 * (r(1) + 35))/(exp(-(r(1) + 35)/10) - 1);
    end
    beta_m = 4 * exp(-(r(1) + 60)/18);


    if r(1) == -50 
        alpha_n = 0.1;
    else
        alpha_n = (-0.01 * (r(1) + 50))/(exp(-(r(1) + 50)/10) - 1);
    end
    beta_n = 0.125 * exp(-(r(1) + 60)/80);

    alpha_h = 0.07 * exp(-(r(1) + 60)/20);
    beta_h = 1/(1 + exp(-(r(1)+30)/10));
    
    if t >10 && t < 30
        iext = -3;
    else
        iext = 0;
    end

    result = zeros(4,1); % v,m,h,n
    result(1) = (1/c) * (iext - (g_k_bar * r(4)^4 * (r(1) - e_k)) - (g_na_bar * r(2)^3 * r(3) * (r(1) - e_na)) - (g_l * (r(1) - e_l)) );
    result(2) = (alpha_m * (1 - r(2))) - (beta_m * r(2));   
    result(3) = (alpha_h * (1 - r(3))) - (beta_h * r(3));
    result(4) = (alpha_n * (1 - r(4))) - (beta_n * r(4));
end