%
%  Elliptic Burst model based on Morris-Lecar cell with subcritical Hopf
%                Victor Matveev, November 23 2015
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

par = [0.04 4.4 2 30 110 0.04 0.0167 1.00];   % Hopf

phi  = par(1);
gCa  = par(2);
V3   = par(3);
V4   = par(4);

Iapp = par(5);
Eps  = par(6)*0.04;
mu   = par(7)*0.3;
gKCa = par(8)*1.2;    

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

ICa    = @(V)      gCa*minf(V).*(V-VCa);
IKCa   = @(V,Ca)   gKCa*Ca./(Ca+1).*(V-VK);
Itotal = @(V,n,Ca) Iapp - ICa(V) - gK*n.*(V-VK) - gL*(V-VL) - IKCa(V, Ca);
f      = @(t,y)    [ Itotal(y(1),y(2),y(3))/C; r(y(1)).*(ninf(y(1))-y(2)); Eps*(-mu*ICa(y(1))-y(3)) ];

options = odeset('RelTol',1e-6);
V0 = VL;

IC = [3 0.1 0.3];
[T Y] = ode45(f, [0, 2*Tmax], IC, options);
V = Y(:,1);

subplot(3, 2, 2); 
plot(T, V, 'r-', 'linewidth', lw); 
ylabel('V(t)');  title('Potential, V(t)');
axis tight;

subplot(3, 2, 4); 
plot(T, Y(:,3), 'r-', 'linewidth', lw);
ylabel('Ca(t)'); title('Calcium concentration, [Ca]');
axis tight;

subplot(3, 2, 6); 
plot(T, IKCa(V,Y(:,3)), 'r-', 'linewidth', lw); 
xlabel('time'); title('I_{K,Ca}(t)');
axis tight;

subplot(3, 2, [1 3 5]); hold on;
title('Elliptic Burst');
n = round(3*length(T)/4);
plot3(Y(n:end,3), Y(n:end,2), Y(n:end,1)); 
xlabel('Ca'); ylabel('n'); zlabel('V');


Vss = []; Nss = []; Vss = []; Css = [];
ic = -25; 
for Ca = linspace(0.5, 0.25, 80)
    h = @(v) Itotal(v, ninf(v), Ca);
    vss = fsolve( h, ic );
    ic = vss;
    Vss = [Vss, vss]; Nss = [Nss, ninf(vss) ]; Css = [Css, Ca];
end
plot3(Css, Nss, Vss, 'r-');







    
