%
%  Morris-Lecar model of an excitable cell (barnacle muscle)
%             Victor Matveev, June 10 2010
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
figure('Position',[1 200 900 450]);  % position the window
lw = 0.5;
Tmax = 200;

V1  = -1.2; V2 = 18;
gK  = 8;    gL = 2;
VCa = 120;  VK = -84; VL=-60;
C = 20;

%par = [0.04 4.4 2 30];   % Hopf
%par = [0.067 4 12 17.4]; % SNIC
par = [0.23 4 12 17.4];   % HOB 

phi = par(1);
gCa = par(2);
V3  = par(3);
V4  = par(4);
    
for Iapp = 33 : 0.1 : 40
    
    minf = @(V) (1 + tanh((V-V1)/V2) )/2;
    ninf = @(V) (1 + tanh((V-V3)/V4) )/2;
    r    = @(V) phi*cosh((V-V3)/(2*V4));

    Itotal = @(V,n) Iapp-gCa*minf(V).*(V-VCa)-gK*n.*(V-VK)-gL*(V-VL);
    f = @(t,y) [ Itotal(y(1),y(2))/C; r(y(1)).*(ninf(y(1))-y(2)) ];

    options = odeset('RelTol',1e-6);
    V0 = VL;
    
    IC = [3 0.1];
    [T Y] = ode45(f, [0, 2*Tmax], IC, options);
    V = Y(:,1);

    i = 1;
    vRange = -75:5:75;
    n = (Iapp-gCa*minf(vRange).*(vRange-VCa)-gL*(vRange-VL)) ./ (gK * (vRange-VK));
 
    subplot(3, 1, [2, 3]); % Phase plane plot goes here
    hold off;
    plot(Y(:,1), Y(:,2),       'r-', 'linewidth', lw); hold on;
    plot(vRange, n,            'k-', 'linewidth', lw); 
    plot(vRange, ninf(vRange), 'b-', 'linewidth', lw);   

    xlabel('V(t)'); ylabel('n(t)'); title('Phase space');
    axis([-75 75 -0.1 0.45]);
    vv = linspace(-75,   75, 10);
    nn = linspace(-0.1, 0.45, 10);
    [VV, NN] = meshgrid(vv, nn);
    DV = (Itotal(VV, NN))/C;
    DN = r(VV).*(ninf(VV)-NN);
    quiver(VV, NN, DV, DN, 'ShowArrowHead', 'off');
    legend('Trajectory', 'V-nullcline', 'w-nullcline', 'Flow field', 'Location', 'SE');

    subplot(3, 1, 1); hold off; % V(t) plot
    plot(T, V, 'r-', 'linewidth', lw); hold on;
    xlabel('time'); ylabel('V(t)');
    axis tight;
    title(['I=', num2str(Iapp)]);

    drawnow;
end
    
