%   Generates Figure 5 
%
%   Oh & Matveev (2010)

clear;
global C ECa EK EL Esyn gCa gK gL Iapp tSyn tRise gSyn Vth TT VV;

gSyn = 1.5; 
fname = ['PRCp4_gSyn',num2str(gSyn),'_tSyn1.mat'];
load(fname);

Tplot = 70;

C = 2; Iapp = 14;
ECa = 120; EK = -84; EL = -60; Esyn = -80;
gCa = 4; gK = 8; gL = 2;
tSyn = 1; tRise = 0.2; Vth = -3;

minf = @(V) (1 + tanh((V+12)/18) )/2; winf = @(V) (1 + tanh((V+8)/6) )/2;
r    = @(V) 0.666*cosh((V+8)/12);     sig  = @(V) (1 + tanh(4*V)) /2;

fpre = @(t,y) [ (-gCa*minf(y(1)).*(y(1)-ECa)-gK*y(2).*(y(1)-EK)-gL*(y(1)-EL)+Iapp)/C; ...
                 r(y(1)).*(winf(y(1))-y(2)); ...
                 -y(3)/tSyn * sig(Vth - y(1)) + (1-y(3))/tRise * sig(y(1) - Vth); ];

options = odeset('RelTol',1e-8);          
[T Y] = ode15s(fpre, [0 400], [-60 0 0]);
[T Y] = ode15s(fpre, [0 200], Y(end,:), options );

[Tmx, Vmx, Imx, Wmx] = getMinMaxP4(T, Y,  1);
[Tmn, Vmn, Imn, Wmn] = getMinMaxP4(T, Y, -1);
P = diff(Tmx);
T0 = P(end);

imin = Imn(end - 2);
ds = diff(Y(imin:end,3));
incr = find(ds > 0, 1) + imin;
iStart = find(Y(incr:end,3) > 1e-7, 1) + incr;
iEnd = find(Y(iStart + 50 : end, 3) < 1e-8, 1) + iStart + 50;
TT = T(iStart:iEnd) - T(iStart); VV = Y(iStart:iEnd, 1);
TT = [TT; TT(end)+1; 200]; VV = [VV; VV(end); VV(end)-10];
PP = spline(TT, VV);

gS = @(t) PPVAL(PP,t) ;

fpost = @(t,y) [ (-gCa*minf(y(1)).*(y(1)-ECa)-gK*y(2).*(y(1)-EK)-gL*(y(1)-EL)-gSyn*y(3).*(y(1)-Esyn)+Iapp)/C; ...
                  r(y(1)).*(winf(y(1))-y(2)); ...
                  -y(3)/tSyn * sig(Vth - gS(t)) + (1-y(3))/tRise * sig(gS(t) - Vth);];

iStart = iStart - 100;        % Just to make a nicer figure with whole AP

Tmin = Tmn( find(Tmn > T(iStart), 1) );
Tmax = Tmx( find(Tmx  > T(iStart), 1) );
              
DT0 = Tmin - Tmax             
DT1 = 0.3; DT2 = 4; DT3 = 5;

tnew = Tmin + DT1;            % First pulse DT1 after V minimum
IC1 = spline(T,Y(:,1),tnew);
IC2 = spline(T,Y(:,2),tnew);
i1  = round(spline(T,1:length(T),tnew));
Ttot = squeeze(T(iStart:i1)) - T(iStart); Vtot = squeeze(Y(iStart:i1,1));

[X Z] = ode15s(fpost, [0 DT2], [IC1 IC2 0], options );
Ttot = cat(1,Ttot,X + Ttot(end)); Vtot = cat(1,Vtot,Z(:,1));
IC = Z(end,1:3);
[X Z] = ode15s(fpost, [0 DT3], IC, options );
Ttot = cat(1,Ttot,X + Ttot(end)); Vtot = cat(1,Vtot,Z(:,1));
IC = Z(end,1:3);
[X Z] = ode15s(fpost, [0 1.8*T0], IC, options );
Ttot = cat(1,Ttot,X + Ttot(end)); Vtot = cat(1,Vtot,Z(:,1));
[T2 Y2] = ode15s(fpre, [0 1.9*T0], Y(iStart,:), options );

subplot(211);
plot(Ttot,Vtot,'k-','linewidth',2); hold on;
plot(T2,Y2(:,1),'k-','linewidth',1);
ylabel('V(t), mV');
set(gca,'box','off');
axis([0 Tplot -60 20]);

[Tm, Vm, I, Wm] = getMinMaxP4(X, Z,  1);
Actual = ( DT0 + DT1 + DT2 + DT3 + Tm(2) ) / T0 - 1

phi1 = (DT0 + DT1)/T0;
Delta1 = spline(phi, PRC, phi1);
phi2 = phi1 - Delta1 + DT2/T0;

if phi2 >= 0
    Delta2 = spline(phi, PRC, phi2);
else
    Delta2 = spline(phiN, PRCN, phi2);
end

phi3 = phi2 - Delta2 + DT3/T0;

if phi3 >= 0
    Delta3 = spline(phi, PRC, phi3);
else
    Delta3 = spline(phiN, PRCN, phi3);
end

Predicted = Delta1 + Delta2 + Delta3

t0 = Tmax - T(iStart); ph0=1-t0/T0;
t1 = t0+DT0+DT1;
t2 = t0+DT0+DT1+DT2;
t3 = t0+DT0+DT1+DT2+DT3;

tt = [0, t0, t0, t1, t1, t2, t2, t3, t3, t3+Tm(2), t3+Tm(2), Tplot];
phase = [ph0, 1, 0, phi1, phi1-Delta1, phi2, phi2-Delta2, phi3, phi3-Delta3, 1, 0, (Tplot-(t3+Tm(2)))/T0]; 

tt0 = [t0, t0+T0, t0+T0, t0+2*T0, t0+2*T0];
pp0 = [0, 1, 0, 1, 0]; 

[Tm, Vm, I, Wm] = getMinMaxP4(T2, Y2, 1);
arrow([t0 0], [Tm(end) 0], 'Ends', 'both', 'Length', 12);
text(t0 + Tm(end)/2, -5, 'T_0');
arrow([t0 10], [tt(end-1) 10], 'Ends', 'both', 'Length', 12);
text((t0+tt(end-1))/2,  15, 'T_3');

subplot(212);
pfs = 13;
plot(tt,phase,'k-','linewidth',2); hold on;
plot(tt0,pp0,'k:','linewidth',1);
hold on;
axis([0 Tplot min(phase)-0.1 1.2]);
plot([0 Tplot], [1 1], 'k:');
plot([0 Tplot], [0 0], 'k:');

plot(t1, phi1, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);
plot(t2, phi2, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);
plot(t3, phi3, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);

h=text(t1+1, phi1+0.1, '\phi_1', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_1$');
h=text(t2+1, phi2+0.1, '\phi_2', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_2$');
h=text(t3+1, phi3+0.1, '\phi_3', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_3$');

plot([t1 t1],[min(phase) 0.75],':k');
plot([t2 t2],[min(phase) 0.75],':k');
plot([t3 t3],[min(phase) 0.75],':k');

arrow([t1 0.7], [t2 0.7], 'Ends', 'both', 'Length', 9);
arrow([t2 0.7], [t3 0.7], 'Ends', 'both', 'Length', 9);
text(0.5*(t1 + t2)-1, 0.85, '\Deltat_1');
text(0.5*(t2 + t3)-1, 0.85, '\Deltat_2');

xlabel('time, ms'); 
h = ylabel('\phi(t)','fontsize',pfs);
%set(h,'Interpreter','latex','string','$\varphi(t)$');
