% ### EXhhFoundation.m ###      2018.03.19      C. Bergevin

% [IN PROGRESS; see "To DO" below]
% GOAL: Illustrate rationale behind Hodgkin-Huxley (1952e) bits (re their Fig.2):
% O "If g_k is used as a variable the end of the record can be fitted by a first-order 
% equation but a third- or fourth-order equation is needed to describe the
% beginning"
% o "A useful simplification is achieved by supposing that g_k is
% proportional to the fourth power of a variable which obeys a first-order equation. 
% [...] The rise in conductancetherefore shows a marked inflexion, while the fall is a simple
% exponential." 
% o [from Nelson] 
% "Hodgkin and Huxley realized that a better fit could be obtained if they considered the 
% conductance to be proportional to a higher power of n. Figure 3 shows the results of fitting the 
% conductance data using a form [eqn] with powers of j ranging from 1 to 4. Using this sort of 
% fitting procedure, Hodgkin and Huxley determined that a reasonable fit to the K+ conductance data 
% could be obtained using an exponent of j=4. Thus they arrived at a description for the K+ 
% conductance under voltage clamp conditions"

% --- Equations
% o Consider the generic 1st order ODE from H&H 1952e (eqns 7-10)
% [using x as a dummy variable name for n, m, or h]
% + \dot{x} = \alpha_x*(1-x) - beta_x*n
% OR AS
% + \tau*\dot{x} + x = x_{infty}
% o where \tau = 1/(\alpha_x + \beta_x) and x_{infty} = \alpha_x/(\alpha_x + \beta_x)
% or \alpha_x= x_{infty}/\tau and \beta_x= (1-x_{infty})/\tau

% o Solutions are of form:
% + x= x_{infty}+ (x_o - x_{infty})* exp(-t/\tau)
% where x_o is the initial condition, i.e., x_o= x(t=0)

% --- To DO
% o something seems fishy re Weiss: I can't get a decent fit to HH data 
% (g_k for Vm= 88 mV, see Weiss 1996 vol.2 Fig.4.23) using 

% o do actual nonlinear regression for different models to those data and
% determine best fit params.
% --> better motivates Fig.1 for comparing n^j for various j
% --> possible to determine appropriate values of n_infty and tauN?? (and
% compare to analytic expressions)
% o repeat both items above for the different membrane polarizations and 
% determine the functional form of the how they depend upon Vm 
% --> presumably should match/reproduce Weiss 1996 vol.2 Fig.4.25 and
% eqns.4.14-4.20

% --- REFs
% o H&H 1952e = J. Physiol. (1952) 117, 500-544
% o Nelson, M.E. (2004) "Electrophysiological Models In: Databasing the 
% Brain: From Data to Knowledge" 

clear
% ==============================================================

% --- generic variable "x" params. (Fig.1)
tau= 4;  % [ms]
xINFTY= 1;
x0= 0;      % initial condition

% --- K (potassium) params.
if (1==0)  % Weiss eqns.4.18 and 4.19 (pg.191)
    GkBAR= 36;
    aK= @(X) -0.01*(X+50)./(exp(-0.1*(X+50))-1);
    bK= @(X) 0.125* exp(-0.0125*(X+60));
    Vm= 28;  % membrane potential for Gk vals below
    Vmx= linspace(-100,75,500);  % functional form (Weiss uses "membrane potential")
else  % H&H 1952 eqns. 12 & 13 (and table 1 caption)
    GkBAR= 24.31; % val. listed (as "figure") in Table 1 caption(?)
    %GkBAR= 21.0;  % better fit value!
    aK= @(X) 0.01*(X+10)./(exp((X+10)/10)-1);
    bK= @(X) 0.125* exp(X/80);
    Vm= -88;  % membrane potential for Gk vals below
    Vmx= -linspace(-40,135,500);  % func. form (HH use disp. re rest. pot.)
end
% Gk conductance values extracted (via deplot) for Vf-Vi= 88 mV
% (presumably Vi=-60 mV such that Vm=28 mV)
% (from Weiss 1996 vol.2 Fig.4.23)
Gk= [0.1358   0.0000 
    0.3387    0.2182
    0.5838    1.0731
    0.7695    2.0969
    1.1035    4.3834
    1.5119    7.0410
    2.0068    9.3710
    2.9435   12.5883
    4.2046   15.5825
    6.4358   16.7832
    8.8070   16.8921];

% ==============================================================
% ---
t= linspace(0,20,1000);  % [ms]
% --- Fig.1: Simple visual comparison of powers of 1st order kinetics
x= xINFTY+ (x0 - xINFTY)* exp(-t/tau);
figure(1); clf;
h1= plot(t,x,'LineWidth',2); hold on; grid on;
h2= plot(t,x.^2,'LineWidth',2); h3= plot(t,x.^3,'LineWidth',2); 
h4= plot(t,x.^4,'LineWidth',2); xlabel('Time [ms]'); ylabel('x^j');
legend([h1 h2 h3 h4],'x','x^2','x^3','x^4','Location','SouthEast');
title('First order kinetics (x) to different powers') 
% --- Fig.2: Comparison to actual HH data
tauN= 1/(aK(Vm)+bK(Vm)); % [ms]
nINFTY= aK(Vm)/(aK(Vm)+bK(Vm));
%n= nINFTY*(1- exp(-t/tauN));
n0=0.24;
n= nINFTY+ (n0 - nINFTY)* exp(-t/tauN);
GkMOD= GkBAR* n.^4;
figure(2); clf;
k1= plot(Gk(:,1),Gk(:,2),'ko','MarkerSize',6,'LineWidth',2); hold on; grid on;
k2= plot(t,GkMOD,'r--','LineWidth',2); xlim([0 1.1*max(Gk(:,1))])
xlabel('Time [ms]'); ylabel('G_k [mS/cm^2]');
title('G_k conductance data (Vm= 88 mV) from Weiss vol.2 Fig.4.23 (from Hodgkin 1964, Fig.29)')
% --- Fig.3: Potassium params. re their voltage depend.
alpha= aK(Vmx); beta= bK(Vmx); tauNf= 1./(alpha+beta); nINFTYf= alpha./(alpha+beta);
figure(3); clf;
% NOTE: what precisely the x-axis represents here depends upon whether you
% use the HH convention or that of Weiss
subplot(221); plot(Vmx,alpha); hold on; grid on; xlabel('pot. diff [V]'); ylabel('\alpha_n')
subplot(222); plot(Vmx,beta); hold on; grid on; xlabel('pot. diff [V]'); ylabel('\beta_n')
subplot(223); plot(Vmx,tauNf); hold on; grid on; xlabel('pot. diff [V]'); ylabel('\tau_n')
subplot(224); plot(Vmx,nINFTYf); hold on; grid on; xlabel('pot. diff [V]'); ylabel('n_{\infty}')