% HodgHux.m
%              19 September 2005
% Hodgkin-Huxley simple simulator
clear

% paramaters  for the membrane patch of pi 30^2 um^2 

C = 1 ; % uF/cm^2  membrane capacitance
gg = [36 120] ; % mS/cm^2  K and Na max conductances
gl = 0.3 ; % mS/cm^2  the leak conductance
E = [-12 115 0]' ; % mV , K, Na, leak to be calculated
Vr = -70; % mV resting potential

% Resting conductances (for V = 0)

alf0 = [0.1/(exp(1)-1), 2.5/(exp(2.5)-1), 0.07] ;
bet0 = [0.125 , 4 , 1/(exp(3)+1)] ;

tau0 = 1./(alf0+bet0) ;
nmh0 = alf0.*tau0 ;

% resting conductances: K Na 
g0 = gg.*[nmh0(1)^4 , nmh0(3)*nmh0(2)^3]  ;

% the leak battery
E(3) = -g0*E(1:2)/gl;  %  10.5989

tm = 20;  % simulation time
% time step for integration
dt = 0.01 ;   % 10 us
t = 0:dt:tm ; % time variable
kk = max(size(t)) ;

% injected current Inj = 10 nA
Inj = 10;  % na from t1:t2
t1 = 3;  t2 = t1+100*dt;     % 3:4 ms
Iin = zeros(kk,1) ; 
Iin(round(t1/dt):round(t2/dt)) = Iin(round(t1/dt):round(t2/dt))+Inj ;

% Initial conditions
V = zeros(kk,1) ;       % membrane voltage
x = zeros(kk, 3) ; x(1,:) = nmh0; % activation particles
gKNa = zeros(kk, 2) ; % potassium and sodium conductances

for k = 1:kk-1

 % alpha and beta functions
 Vk = V(k) ;
 alfa = [0.01*(10-Vk)./(exp((10-Vk)/10)-1) % n
         0.1*(25-Vk)./(exp((25-Vk)/10)-1) % m 
         0.07*exp(-Vk/20) ]' ;  % h

 beta = [0.125*exp(-Vk/80) % n
         4*exp(-Vk/18)    % m 
        1./(exp((30-Vk)/10)+1)]' ; % h
    
 % time constants and steady-state values

 tau = 1./(alfa + beta) ;
 xs = alfa.*tau ;

 % conducances
 
 gKNa(k,:) = gg.*[x(k,1)^4 x(k,3)*x(k,2)^3] ;
 g = [gKNa(k,:) gl];

 % activation particles
 x(k+1,:) = x(k,:) + dt*(xs - x(k,:))./tau ;

 % membrane voltage
 V(k+1) = V(k) + (g*(E-V(k)) + Iin(k))*(dt/C) ;

end

 myfigA(1,9,7)  % figure(1)
 plot(t,x), grid on
 xlabel('time [msec]'), ylabel('activation particles')
 text(0.5, x(1,1)+0.03, 'n (K)') 
 text(0.5, x(1,2)+0.03, 'm (Na)')
 text(0.5, x(1,3)+0.03, 'h (Na)')
 title('Activation particles during action potential')
 
  myfigA(2,9,7)  % figure(2)
 plot(t, gKNa), grid on
 xlabel('time [msec]'), ylabel('conductances [mS/cm^2]')
 text(9, 8.1, 'g_K') 
 text(7, 21, 'g_{Na}')
 title('Potassium and Sodium conductances')
 
 
 myfigA(3,9,7)  % figure(3)
 plot(t, [V+Vr Iin-95], '-', [0 20], [1 1]*Vr, '--'), grid on
 text(1.8,-80, 'I_{in}=10nA')
 text(7, 7, 'V_m')
 xlabel('time [msec]'), ylabel('Action potential [mV]')
 title('Injected current and action potential')