function TsynPRC(g, tS)

% Generate accurate PRC (have some patience - this runs slowly)
% Arguments are: "g" - inhibitory synaptic conductance
%                "tS" - synaptic time constant
%    Oh & Matveev (2010)
%

global C ECa EK EL Esyn gCa gK gL Iapp tSyn tRise Vth TT VV PP gSyn;
gSyn = g;
tSyn = tS;

C = 2; Iapp = 14;
ECa = 120; EK = -84; EL = -60; Esyn = -80;
gCa = 4; gK = 8; gL = 2;
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;

% ODE system for the "presynaptic" cell: 

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-10);
          
[T Y] = ode15s(fpre, [0 210], [-60 0 0], options);  % Pre-equilibration run
IC = Y(length(T),:);
[T Y] = ode15s(fpre, [0 1000], IC, options );

[Tm, Vm, I, Wm] = getMinMaxP4(T, Y,  1);           % Extract accurate spike times
P = diff(Tm);
T0 = mean(P(end-4:end));                           % Period = average of last 4 ISIs

% This block defines the limit cycle as the interval between the last two spikes:

IN = I(end-1) - 1;            % First index of limit cycle
IL = I(end)   + 1;            % Last index of limit cycle
Tstart = Tm(end-1);           % Initial time of limit cycle
LCT = T(IN:IL) - Tstart;      % Time array for a single limit cycle
LCY = Y(IN:IL,:);             % Limit cycle trajectory
LCV = spline(LCT, LCY(:,1));  % Voltage trace of the limit cycle
LCW = spline(LCT, LCY(:,2));  % Recovery variable during the limit cycle

[Tm, Vm, Im, Wm] = getMinMaxP4(T, Y, -1);  % Extract accurate times of voltage minimum
Wmin = Wm(end);
Vmin = Vm(end);                            % Value of voltage minimum

imin = Im(end - 2);          % This code block extract voltage waveform of a single spike
ds = diff(Y(imin:end,3));
incr = find(ds > 0, 1) + imin;
iStart = find(Y(incr:end,3) > 1e-7, 1) + incr;                % Find location of spike start
iEnd = find(Y(iStart + 50 : end, 3) < 1e-8, 1) + iStart + 50; % Find location of spike end
TT = T(iStart:iEnd) - T(iStart); VV = Y(iStart:iEnd, 1);      % Spike time array
TT = [TT; TT(end)+1; 200]; VV = [VV; VV(end); VV(end)-10];    % Append a value so that voltage goes below threshold
PP = spline(TT, VV);                                          % Spline for V(t) of a single spike

[tm, vm, im, wm] = getMinMaxP4( T(iStart:iEnd), Y(iStart:iEnd,:), 1 );  % Find phase lag between s(t) onset and max V(t)
phi0 = ( tm(1) - T(iStart) ) / T0                                       % Phase lag = phi0


gS = @(t) PPVAL(PP,t) ;  % PP presynaptic voltage trace is injected into the postsynaptic cell equations below:

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);];

j = 1; tnew = 0; told = tnew; 
step = 0.01 * T0;
minStep = 0.00025 * T0;
maxStep = 0.01 * T0;

% Main loop: calculate PRC as a function of phase phi (predictor-corrector stepping algorithm):

while tnew < T0
    phi(j) = tnew / T0;
    IC1 = PPVAL(LCV, tnew);
    IC2 = PPVAL(LCW, tnew);
    [X Z] = ode15s(fpost, [0 T0*2.5], [IC1 IC2 0], options );
    [Tm, Vm] = getMinMaxP4(X, Z, 1);
        
    for pos = 1 : 2   
        PRC(j) = phi(j) + Tm(pos)/T0 - 1;           % Value of 1-st order PRC
        PRC2(j) = (Tm(pos+1) - Tm(pos)) / T0 - 1;   % Value of 2-nd order PRC (quite small)
        if (Vm(pos) > Vth) 
            break;
        end
    end
    
    if j>2 && abs((PRC(j)-PRC(j-1))/(PRC(j)+PRC(j-1))) > 0.01*(1-log10(PRC(j))) && step > minStep
        step = step * 0.5;
        tnew = told + step;
        fprintf('Correction: time=%g %g\n', PRC(j), PRC(j-1));
        continue;
    else
        fprintf('[%g %g %g] %g %g %g %d\n', told, tnew, step, phi(j), PRC(j), PRC2(j), pos);
        told = tnew;
        tnew = told + step; 
        if step < maxStep
            step = step * 1.1;
        end
        j = j + 1;
    end
end;


% Calculate PRC on the negative phase interval (PRCN as a function of phiN):

j = 1; V = Vmin;

while j == 1 || PRCN(j-1) > 3e-3
    [X Z] = ode15s(fpre, [0 T0*1.5], [V, Wmin, 0], options );
    Tm = getMinMaxP4(X, Z, 1);
    if Tm(1) < 0.9 * T0
        break
    end;
    phiN(j) = 1 - Tm(1)/T0;
    
    [X Z] = ode15s(fpost, [0 T0*1.5], [V, Wmin, 0], options );
    [Tm, Vm] = getMinMaxP4(X, Z, 1);
    
    for pos = 1 : 2
        PRCN(j) = phiN(j) + Tm(pos)/T0 - 1;
        if (Vm(pos) > Vth) 
            break;
        end
    end
    fprintf('%g %g\n', phiN(j), PRCN(j));
    j = j + 1;
    V = V - j/10;
end;

fname = ['PRCp4_gSyn',num2str(gSyn),'_tSyn',num2str(tSyn),'.mat'];
save(fname,'phi','PRC','PRC2','phiN','PRCN','phi0','LCT','LCV','PP','LCW','T0');


