function hw3p2
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology
%
% hw3p2.m - Homework #3, Problem #2
%
% HODGKIN-HUXLEY: Response of Squid Axon to Current Injection
%
% M. Nelson  1995, 1997
%----------------------------------------------------------------------

clear
clf
orient landscape

%%%%%%%%%%%%%%%%%%%%%%%
% INITIALIZATION
%%%%%%%%%%%%%%%%%%%%%%%
DT = 0.025;
TMAX = 50;
t = 0:DT:TMAX;

VREST = 0.0;

GNA = 120.0;	% mS/cm^2
GK = 36.0;	% mS/cm^2
GLEAK = 0.3;	% mS/cm^2

ENA = 115.0;	% mV
EK = -12.0;	% mV
ELEAK = 10.613;	% mV

C_M = 1.0;	% uF/cm^2

INJECT_LIST = [1.0 10.0 -1.0 -10.0];

v = zeros(size(t));
m = zeros(size(t));
n = zeros(size(t));
h = zeros(size(t));

for ival = 1:4

    I = INJECT_LIST(ival)*(t>5 & t<=30);

    %%%%%%%%%%%%%%%%%%%%%%%
    % MAIN SIMULATION LOOP
    %%%%%%%%%%%%%%%%%%%%%%%

    v(1) = VREST;						% initial v
    m(1) = alpha_m(VREST)/(alpha_m(VREST)+beta_m(VREST));	% initial m
    h(1) = alpha_h(VREST)/(alpha_h(VREST)+beta_h(VREST));	% initial h
    n(1) = alpha_n(VREST)/(alpha_n(VREST)+beta_n(VREST));	% initial n

    for i=2:length(t)

	M = m(i-1);	% get values from last time step
	H = h(i-1);	% (hopefully this substitution makes 
	N = n(i-1);	% the following code a bit easier to read)
	V = v(i-1);

	gNa = GNA * M^3 * H;
	gK  = GK  * N^4;

	mdot = alpha_m(V)*(1-M) - beta_m(V)*M;
	hdot = alpha_h(V)*(1-H) - beta_h(V)*H;
	ndot = alpha_n(V)*(1-N) - beta_n(V)*N;
	vdot = (I(i-1) - gNa*(V-ENA) - gK*(V-EK) - GLEAK*(V-ELEAK))/C_M;

	m(i) = m(i-1) + mdot*DT;		% Euler integration
	h(i) = h(i-1) + hdot*DT;
	n(i) = n(i-1) + ndot*DT;
	v(i) = v(i-1) + vdot*DT;
    end

    subplot(2,2,ival);
    plot(t,v,'b-'); hold on;
    plot(t,I,'r:')
    str = sprintf('Inject = %2.2f \\muA/cm^2',INJECT_LIST(ival));
    title(str);
    drawnow;
end

%%%%%%%%%%%%%%%%%%
% label plots
%%%%%%%%%%%%%%%%%%
subplot(221); 
ylabel('v (mV)')
axis([0 TMAX -2 2])

subplot(223); 
xlabel('t (msec)'); ylabel('v (mV)')
axis([0 TMAX -2 2])

subplot(222); 
ylabel('v (mV)')
axis([0 TMAX -50 150])

subplot(224); 
xlabel('t (msec)'); ylabel('v (mV)')
axis([0 TMAX -50 150])

%%%%%%%%%%%
% alpha_h
%%%%%%%%%%%
function rate = alpha_h(v)
  rate = 0.07*exp(-v/20.0);

%%%%%%%%%%%
% alpha_m
%%%%%%%%%%%
function rate = alpha_m(v)
eps = 1e-12;

  rate = zeros(size(v));		% initialization

  idx = find(abs(v-25) > eps);		% HANDLE NORMAL (NON-SINGULAR) CASE
  if(~isempty(idx))
    rate(idx) = 0.1*(-v(idx)+25.0) ./ (exp((-v(idx)+25.0)/10.0)-1.0);
  end

  idx = find(abs(v-25) <= eps);		% HANDLE SIGULARITY AT V=25.0 
  if(~isempty(idx)) 
    rate(idx) = 1.0; 
  end

%%%%%%%%%%%
% alpha_n
%%%%%%%%%%%
function rate = alpha_n(v)
eps = 1e-12;

  rate = zeros(size(v));

  idx = find(abs(v-10) > eps);
  if(~isempty(idx))
    rate(idx) = 0.01*(-v(idx)+10.0) ./ (exp((-v(idx)+10.0)/10.0) - 1.0);
  end

  idx = find(abs(v-10) <= eps);		% HANDLE SIGULARITY AT V=10.0 
  if(~isempty(idx))
    rate(idx) = 0.1;
  end
       
%%%%%%%%%%%
% beta_h
%%%%%%%%%%%
function rate = beta_h(v)
  rate =  1.0 ./ (exp((-v+30.0)/10.0) + 1.0);

%%%%%%%%%%%
% beta_m
%%%%%%%%%%%
function rate =  beta_m(v)
  rate = 4.0*exp(-v/18.0);

%%%%%%%%%%%
% beta_n
%%%%%%%%%%%
function rate = beta_n(v)
  rate = 0.125*exp(-v/80.0);
