%
%  Tri-stability in the Morris-Lecar model near the 
%      Saddle Homoclinic Orbit Bifurcation
%          Victor Matveev, June 21 2010
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

global C V1 V2 V3 V4 gCa gK gL VK VL VCa Iapp vv;

V1  = -1.2; V2 = 18;  V3 = 12; V4 = 17.4;
gCa = 4;  gK = 8;   gL = 2;
VCa = 120;  VK = -84; VL=-60;
C = 20; Iapp = 37.5;

minf = @(V) (1 + tanh((V-V1)/V2) )/2;
winf = @(V) (1 + tanh((V-V3)/V4) )/2;
r    = @(V) cosh((V-V3)/(2*V4)) * 0.23;

Itotal = @(V,w) -gCa*minf(V).*(V-VCa)-gK*w.*(V-VK)-gL*(V-VL)+Iapp;
f = @(t,y) [ Itotal(y(1),y(2))/C; r(y(1)).*(winf(y(1))-y(2)) ];
          
options = odeset('RelTol',1e-6);
T = 1000;
IC1 = [-60 0];  [T1 Y1] = ode15s(f, [0 T], IC1, options);
IC2 = [0 0.1];  [T2 Y2] = ode15s(f, [0 T], IC2, options);
IC3 = [6 0.3];  [T3 Y3] = ode15s(f, [0 T], IC3, options);

subplot(3,2,2);                  % V(t) plot
hold off;
plot(T1, Y1(:,1), 'b-');
title('Trajectory near low-voltage stable fixed point ');

subplot(3,2,4);                   % Plot activation gates
hold off;
plot(T2, Y2(:,1), 'g-');
title('Trajectory close to the stable limit cycle');

subplot(3,2,6);                   % Plot currents
hold off;
plot(T3, Y3(:,1), 'm-');
title('Trajectory near high-voltage stable fixed point');
xlabel('time'); 

subplot(3,2,[1 3 5]); hold off;    % Phase plane plot
plot(Y1(:,1), Y1(:,2), 'b-', 'linewidth', 1.5); hold on;
plot(Y2(:,1), Y2(:,2), 'g-', 'linewidth', 1.5);
plot(Y3(:,1), Y3(:,2), 'm-', 'linewidth', 1.5); 
hold on;

i = 1;
vRange = -75:5:50;

for vv = vRange   % This calculates the V-nullcline
 w(i) = fsolve(@(w) Itotal(vv,w), 0);
 i = i + 1;
end;

plot(vRange, w, '-', 'color', [0.6 0.6 0.6]);           % plot the V-nullcline
plot(vRange, winf(vRange),'-', 'color', [0.6 0.6 0.6]); % Plot the w-nullcline
xlabel('V(t)'); ylabel('w(t)'); 
title(['Phase space, I_{app}=', num2str(Iapp)] );
axis([-75 50 -0.1 0.5]);
legend('Trajectory 1', 'Trajectory 2', 'Trajectory 3', 'Location', 'SE');

drawnow;