% Generates Figure 7 - Spike-time sequence prediction for a three-cell network
% This program uses both 1st and 2nd order PRCs
%
%   Oh & Matveev (2010)


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

figure(1);
gSyn = 0.15; 
tSyn = 1; 
fname = ['PRCp4_gSyn',num2str(gSyn),'_tSyn',num2str(tSyn),'.mat'];
load(fname);
I = find(diff(PRC2) < 0); J = find(PRC2 < 0);
PRC2(1:I(end)) = 0; PRC2(J) = 0;

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;

ICa  = @(V)     -gCa*minf(V).*(V-ECa);
IK   = @(w,V)   -gK*w.*(V-EK);
IL   = @(V)     -gL*(V-EL);
Isyn = @(s,V)   -gSyn*s.*(V-Esyn);

three = @(t,y) [ (ICa(y(1)) + IK(y(2),y(1)) + IL(y(1)) + Isyn(y(6),y(1)) + Isyn(y(9),y(1)) + 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); ...
                 ...
                 (ICa(y(4)) + IK(y(5),y(4)) + IL(y(4)) + Isyn(y(3),y(4)) + Isyn(y(9),y(4)) + Iapp)/C; ...
                 r(y(4)).*(winf(y(4))-y(5)); ...
                 -y(6)/tSyn * sig(Vth - y(4)) + (1-y(6))/tRise * sig(y(4) - Vth); ...
                 ...
                 (ICa(y(7)) + IK(y(8),y(7)) + IL(y(7)) + Isyn(y(3),y(7)) + Isyn(y(6),y(7)) + Iapp)/C; ...
                 r(y(7)).*(winf(y(7))-y(8)); ...
                 -y(9)/tSyn * sig(Vth - y(7)) + (1-y(9))/tRise * sig(y(7) - Vth)  ];


phases = [1/4 1/2 3/4];
%phases = [0.35 1/2 3/4];

tt = T0*phases;
V = PPVAL(LCV,tt);  W = PPVAL(LCW,tt);  

options = odeset('RelTol',1e-12);          
[T Y] = ode15s(three, [0 500], [V(1) W(1) 0 V(2) W(2) 0 V(3) W(3) 0], options);
subplot(3,1,[1 2]);
plot(T,Y(:,1),'b-'); hold on;
plot(T,Y(:,4),'r-');
plot(T,Y(:,7),'k-');

cls = [0 0 1; 1 0 0; 0 0 0];
pSpike = 1 - phi0;
reset2 = [0 0 0];
time = 0;
tSpike = [1 1];

for i = 1 : 25
    [FPT I] = min(pSpike-phases);
    fprintf('%d T=%4g ',i, FPT);
    phases = phases + FPT;
    tSpike = tSpike + FPT;
    time   = time + FPT;
    subplot(3,1,[1 2]);
    %plot(time*T0, PPVAL(LCV, pSpike*T0), '.', 'color', cls(I,:) );
    
    J = find([1 2 3] ~= I); 
    fprintf('[%g %g]',phases(J));
    
    for j = J
        %plot(time*T0, PPVAL(LCV, phases(j)*T0), '.', 'color', cls(j,:) );
        if tSpike(j) < 0
            fprintf('=R=');
            firstOrder(j) = spline(phi,PRC2,1+tSpike(j)) + spline(phi,PRC,1+tSpike(j)) ;
            secondOrder(j) = 0;
        elseif phases(j) > phiN(1) || ( (tSpike(j) < phiN(1)) && phases(j) > 0)
            K = find(J ~= j);
            k = J(K);
            if tSpike(k)*T0 < 3
                fprintf('+D+');
            end
            firstOrder(j)  = spline(phi,PRC,phases(j));
            secondOrder(j) = spline(phi,PRC2,phases(j));
        else
            fprintf('*N*');
            firstOrder(j)  = spline(phiN,PRCN,phases(j));
            secondOrder(j) = 0;
        end
        reset2(j) = reset2(j) + secondOrder(j);
        phases(j) = phases(j) - firstOrder(j);
        fprintf('(%g', firstOrder(j));
        if secondOrder(j) < 1e-6
            secondOrder(j) = 0.0;
            fprintf('}');
        else
            fprintf(' %g}', secondOrder(j));
        end
        
    end
    
    phases(I) = -reset2(I)-phi0;
    reset2(I) = 0;
    tSpike(I) = -phi0;
    
    fprintf('[%g %g %g]\n',phases);
    
    subplot(313);
    plot([time+phi0, time+phi0] * T0, [0 1], '-', 'color', cls(I,:) ); hold on;
end
   
subplot(3,1,[1 2]);
axis([0 T0*(time+phi0+.1) -50 20]);
set(gca,'box','off');
set(gca,'TickLength',[0.003 0.003]);
set(gca,'XTick',[]);
ylabel('V(t)');

subplot(313);
axis([0 T0*(time+phi0+.1) 0 1]);
set(gca,'box','off');
set(gca,'TickLength',[0.003 0.003]);
set(gca,'YTick',[]);
xlabel('Time');


