function hw3p1
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology 
%
% hw3p1.m - Homework #3, Problem #1
%
% Plot steady-state values and time constants for 
% Hodgkin-Huxley state variables: m,h,n
%
% M. Nelson  1995, 1997
%----------------------------------------------------------------------

clear;
clf;

%--------------------------------
% Define voltage range for plots
%--------------------------------

VMIN =  -50;
VMAX = +150;
DV = 5;

v = VMIN:DV:VMAX;

%-------------
%     m
%-------------
a=alpha_m(v); 
b=beta_m(v);
tau = 1 ./ (a+b);
ss = a .* tau;

subplot(231); plot(v,tau); grid on;
axis('square'); axis([VMIN VMAX 0 1])
xlabel('v (mV)'); ylabel('tau (msec)'); title('\tau_m','fontsize',16);

subplot(234); plot(v,ss); grid on;
axis('square'); axis([VMIN VMAX 0 1])
xlabel('v (mV)'); title('m_\infty','fontsize',16);


%-------------
%     h
%-------------

a=alpha_h(v);
b=beta_h(v);
tau = 1 ./ (a+b);
ss = a .* tau;

subplot(232); plot(v,tau); grid on;
axis('square'); axis([VMIN VMAX 0 10])
xlabel('v (mV)'); ylabel('tau (msec)'); title('\tau_h','fontsize',16);

subplot(235); plot(v,ss); grid on;
axis('square'); axis([VMIN VMAX 0 1])
xlabel('v (mV)'); title('h_\infty','fontsize',16);

%-------------
%     n
%-------------

a=alpha_n(v);
b=beta_n(v);
tau = 1 ./ (a+b);
ss = a .* tau;

subplot(233); plot(v,tau); grid on;
axis('square'); axis([VMIN VMAX 0 10])
xlabel('v (mV)'); ylabel('tau (msec)'); title('\tau_n','fontsize',16);

subplot(236); plot(v,ss); grid on;
axis('square'); axis([VMIN VMAX 0 1])
xlabel('v (mV)'); title('n_\infty','fontsize',16);

orient landscape

%%%%%%%%%%%
% 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-1) > 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);
