%
%  2D reduction of Hodgkin-Huxley Model
%     Victor Matveev, November 2015
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
gNa  = 120; gK = 36;  gL = 0.3;
VNa  = 55;  VK = -77; VL = -54.4;
C    = 5;   I0 = 0;
Iapp = 6;   k  = 0.1;

IC   = [-65.1223415; 0.32];  % Initial condition: [ V(0), n(0) ]

an = @(V) 0.01 * (-(55 + V)) ./ (exp(-0.1*(55 + V)) - 1);
bn = @(V) 0.125 * exp(-(65+V)/80);

am = @(V) -0.1 * (40 + V) ./ (exp(-0.1*(40 + V)) - 1);
bm = @(V) 4 * exp(-(65+V)/18);

minf = @(V) am(V) ./ (am(V) + bm(V));
ninf = @(V) an(V) ./ (an(V) + bn(V));

% Here is the 2D reduction of the Hodgkin-Huxley model:
Itotal = @(V,n) Iapp - gL*(V-VL) - gK*n.^4.*(V - VK) ...
                     - gNa*minf(V).^3.*(0.89 - 1.1*n).*(V-VNa);

f = @(t,y) [ Itotal(y(1), y(2)) / C; an(y(1)).*(1-y(2))-bn(y(1)).*y(2) ];
           
[T, Y] = ode15s(f, [0 60], IC, odeset('RelTol', 1e-6));


subplot(311);
plot(T,Y(:,1));
xlabel('t'); ylabel('V(t)');
title('Membrane potential vs time');

subplot(3,1,[2 3]);
plot(Y(:,1), Y(:,2)); hold on;

VV = -76 : sqrt(1.1) : 50;  % V-values for the nullclines and phase plane
nn = zeros(size(VV));
ic = 0;

for k = 1 : length(VV)  % Build the V-nullcline
   V = VV(k);
   g = @(n) Itotal(V, n);
   nn(k) = fsolve( g, ic);
   ic = nn(k);
end

plot(VV, nn,       'r-');
plot(VV, ninf(VV), 'b-');

xlabel('t'); %ylabel('V(t)');
title('Gating variables vs time');
      
vv = linspace(-75,  55, 10);
nn = linspace(0.2, 0.8, 10);
[VV, NN] = meshgrid(vv, nn);
DV = Itotal(VV,NN) / C;
DN = an(VV).*(1-NN)-NN.*bn(VV);

quiver(VV, NN, DV, DN, 'ShowArrowHead', 'off');
axis([-80 60 0.2 0.8]);

