
% This is an auxiliary script (subroutine) for script Fig6_ThreePulseVarGsyn.m
% This script uses the PRC generated by PRC_Tsyn.m to predict the phase delay produced by three spikes
%    Oh & Matveev (2010)

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

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;

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) ];
           
Tmin = phiN(1)*T0;
DT0  = Tmin; 
tnew = Tmin + DT1;            % First pulse DT1 after V minimum
IC1  = PPVAL(LCV,tnew);
IC2  = PPVAL(LCW,tnew);

options = odeset('RelTol',1e-9);
[X Z]   = ode15s(fpost, [0 DT2], [IC1 IC2 0], options );
IC      = Z(end,1:3);
[X Z]   = ode15s(fpost, [0 DT3], IC, options );
IC      = Z(end,1:3);
[X Z]   = ode15s(fpost, [0 1.8*T0], IC, options );

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

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

if phi2 >= 0
    Delta2 = spline(phi, PRC, phi2);
    Diff2  = spline(phi, PRC, phi2+1e-6)-Delta2;
    Same2  = Delta2; First2 = Delta2; Zero2 = Delta2;
else
    Delta2 = spline(phiN, PRCN, phi2);
    Diff2  = spline(phiN, PRCN, phi2+1e-6)-Delta2;
    Same2  = Delta1;
    First2 = spline(phi, PRC, phi0);
    Zero2  = 0;
end

phi3 = phi2 - Delta2 + DT3/T0;

if phi3 >= 0
    Delta3 = spline(phi, PRC, phi3);
    Diff3  = spline(phi, PRC, phi3+1e-6)-Delta3;
    Same3  = Delta3; First3 = Delta3; Zero3 = Delta3;
else
    Delta3 = spline(phiN, PRCN, phi3);
    Diff3  = spline(phiN, PRCN, phi3+1e-6)-Delta3;
    First3 = spline(phi, PRC, phi0);
    Zero3  = 0;
    if phi2 < 0
     Same3 = Same2;
    else
     Same3 = Delta2;
    end;
end

Predicted = Delta1 + Delta2 + Delta3
Same      = Delta1 + Same2  + Same3
First     = Delta1 + First2 + First3
Zero      = Delta1 + Zero2  + Zero3
Slope     = (Diff1 + Diff2  + Diff3) / 1e-6;
