%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BC2008_fig1
%
% This script runs in MATLAB (and perhaps OCTAVE) and generates
% figure 1 in the paper:
%
% Spike-train spectra and network response functions for
% non-linear integrate-and-fire neurons
% Biological Cybernetics 99: 381-392 (2008)
%
% Please reference this publication in any document that uses the code.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COPYRIGHT: Magnus Richardson (2008).
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% See <http://www.gnu.org/licenses/> for a copy of  
% the GNU General Public License
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BC2008_fig1

DT=3; VT=-53; 
tref1=2; tref2=8; tref=tref1+tref2; % tref1 spike-top tref2 post-spike
tau0=20; Vre=-60; Vth=20; dV=0.01;

%E0B=-45; sig0B=2; VlB=-100;  dVB=0.01;
E0B=-50; sig0B=2; VlB=-100;  dVB=0.01;
E0N=-60; sig0N=6; VlN=-100;  dVN=0.01;

th2deg=360/(2*pi);
K=1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the steady state graphs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

E0Bl=-70; E0Bh=-40; dE0B=0.05; E0Ball=[E0Bl:dE0B:E0Bh]'; nE0B=length(E0Ball);
r0Ball=zeros(nE0B,1); 

for k=1:nE0B
    [V,P0,J0,r0Ball(k)]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig0B,E0Ball(k),dV,VlB);
end

E0Nl=-70; E0Nh=-40; dE0N=0.05; E0Nall=[E0Nl:dE0N:E0Nh]'; nE0N=length(E0Nall);
r0Nall=zeros(nE0N,1); 

for k=1:nE0N
    [V,P0,J0,r0Nall(k)]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig0N,E0Nall(k),dV,VlN);
end

E0Zl=VT-DT; E0Zh=-40; dE0Z=0.05; E0Zall=[E0Zl:dE0Z:E0Zh]'; nE0Z=length(E0Zall);
r0Zall=zeros(nE0Z,1); sig0Z=0.01;

for k=1:nE0Z
    [V,P0,J0,r0Zall(k)]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig0Z,E0Zall(k),dV,VlB);
end
r0Zall(1)=0;

[VB,P0B,J0B,r0B]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig0B,E0B,dV,VlB);
[VN,P0N,J0N,r0N]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig0N,E0N,dV,VlN);

subplot(2,3,1);
plot(E0Ball,K*r0Ball,E0Nall,K*r0Nall,E0Zall,K*r0Zall,'-',E0B,K*r0B,'o',E0N,K*r0N,'o');
xlabel('E0 (mV)'); ylabel('Rate (Hz)');

subplot(4,3,7)
plot(VB,K*J0B,'-',VN,K*J0N,'-')
xlabel('V (mV)'); ylabel('J_0(V)');

subplot(4,3,10)
plot(VB,P0B,'-',VN,P0N,'-')
xlabel('V (mV)'); ylabel('Contin. P(V)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the theory for the FPT, STR and SPECTRUM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% f range. f is in kHz
fl=0.001/1000; fh=1000/1000; nfpos=1000;
fpos=exp(linspace(log(fl),log(fh),nfpos))';
f=[-fpos(nfpos:-1:1);0;fpos]; nf=length(f);
fHz=f*1000; w=2*pi*f; 
dw=([diff(w);w(nf)-w(nf-1)]+[w(2)-w(1);diff(w)])/2;
dw0=dw(find(w==0));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the fourier transforms of the quanties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hfptB=zeros(nf,1);  hstrB=zeros(nf,1); hfspB=zeros(nf,1);
hfptN=zeros(nf,1);  hstrN=zeros(nf,1); hfspN=zeros(nf,1);

for m=1:nf
    [hfptB(m),hstrB(m),hfspB(m)]=ftp_str_fsp(tau0,sig0B,E0B,Vth,Vre,VT,DT,tref,dVB,VlB,w(m),dw0);
    [hfptN(m),hstrN(m),hfspN(m)]=ftp_str_fsp(tau0,sig0N,E0N,Vth,Vre,VT,DT,tref,dVN,VlN,w(m),dw0);
end

[t,fptB]=ift(dw,w,hfptB); fptB=real(fptB);
mfptB=sum(fptB.*t)*(t(2)-t(1)); % this is the mean first passage time
[t,fptN]=ift(dw,w,hfptN); fptN=real(fptN);
mfptN=sum(fptN.*t)*(t(2)-t(1)); % this is the mean first passage time

[t,strB]=ift(dw,w,hstrB); strB=real(strB);
[t,strN]=ift(dw,w,hstrN); strN=real(strN);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% finish with some plots

subplot(4,3,5); plot(t,1000*fptB,'b-');
xlabel('Time (ms)'); ylabel('fpt (Hz)');

subplot(4,3,8); plot(t,1000*strB,'k-',[min(t),max(t)],1000*r0B*[1,1],'k:');
xlabel('Time (ms)'); ylabel('str (Hz)');

subplot(4,3,11); plot(fHz,1000*hfspB);
xlabel('Frequency (Hz)'); ylabel('spectrum (Hz)');
axis([-500 500 0 2*r0B*1000]);

subplot(4,3,6); plot(t,1000*fptN,'b-');
xlabel('Time (ms)'); ylabel('fpt (Hz)');

subplot(4,3,9); plot(t,1000*strN,'k-',[min(t),max(t)],1000*r0N*[1,1],'k:');
xlabel('Time (ms)'); ylabel('str (Hz)');

subplot(4,3,12); plot(fHz,1000*hfspN);
xlabel('Frequency (Hz)'); ylabel('spectrum (Hz)');
axis([-500 500 0 2*r0N*1000]);

drawnow;

if 1
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the sim results for the FPT, STR and SPECTRUM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[tisiB,sisiB,tstrB,sstrB,ffspB,sfspB]=sims(tau0,sig0B,E0B,Vth,Vre,VT,DT,tref);
[tisiN,sisiN,tstrN,sstrN,ffspN,sfspN]=sims(tau0,sig0N,E0N,Vth,Vre,VT,DT,tref);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add sims to plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(4,3,5); hold on;  plot(tisiB,1000*sisiB,'.'); hold off;
subplot(4,3,8); hold on;  plot(tstrB,1000*sstrB,'.'); hold off;
subplot(4,3,11); hold on;  plot(fHz,1000*hfspB,ffspB*1000,1000*sfspB,'.'); hold off

subplot(4,3,6); hold on;  plot(tisiN,1000*sisiN,'.'); hold off;
subplot(4,3,9); hold on;  plot(tstrN,1000*sstrN,'.'); hold off;
subplot(4,3,12); hold on;  plot(fHz,1000*hfspN,ffspN*1000,1000*sfspN,'.'); hold off

drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the example trajectories
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=0.1; T=2000;
[tt,VVB]=sim_example(tau0,sig0B,E0B,Vth,Vre,VT,DT,tref,T,dt);
[tt,VVN]=sim_example(tau0,sig0N,E0N,Vth,Vre,VT,DT,tref,T,dt);

subplot(4,3,2); plot(tt,VVB); xlabel('Time (ms)'); ylabel('V_B (mV)');
subplot(4,3,3); plot(tt,VVN); xlabel('Time (ms)'); ylabel('V_N (mV)');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   steady-state r0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,P0,J0,r0]=EIF_r0_thy(DT,VT,Vth,Vre,tau,tref,sig,E0,dV,Vl);

V=[Vl:dV:Vth]'; n=length(V);

F=(V-E0-DT*exp((V-VT)/DT))/sig^2; 
F1=exp(dV*F); 

pos=(find(F==0));   % this bit removes any inconvenient zeros in the denom
if isempty(pos)
    F2=(F1-1)./F;
else
    F(pos)=1;       % this arbitrary
    F2=(F1-1)./F;  
    F2(pos)=dV;
end

kre=find(V==Vre);

j0=zeros(n,1); j0(n)=1; 
p0=zeros(n,1);         

for k=n:-1:2
j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*F1(k) + tau*j0(k)*F2(k)/sig^2;
end

r0=1/(sum(dV*p0)+tref); 
P0=p0*r0; J0=j0*r0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hfpt,hstr,hfsp]=ftp_str_fsp(tau0,sig,E0,Vth,Vre,VT,DT,tref,dV,Vl,w,dw0);

V=[Vl:dV:Vth]'; n=length(V);

F=(V-E0-DT*exp((V-VT)/DT))/sig^2; 
F1=exp(dV*F); 

pos=(find(F==0));   % this bit removes any inconvenient zeros in the denom
if isempty(pos)
    F2=(F1-1)./F;
else
    F(pos)=1;       % this arbitrary
    F2=(F1-1)./F;  
    F2(pos)=dV;
end

kre=find(V==Vre);

% the ftp is calculated, then the str derived, then the fsp

hp0=zeros(n,1); hj0=zeros(n,1);
hpr=zeros(n,1); hjr=zeros(n,1); hjr(n)=1;

for k=n:-1:2
 
    hj0(k-1)=hj0(k) + dV*i*w*hp0(k) - exp(-i*w*tref)*(k==kre+1);
    hp0(k-1)=hp0(k)*F1(k) + tau0*hj0(k)*F2(k)/sig^2;

    hjr(k-1)=hjr(k) + dV*i*w*hpr(k);
    hpr(k-1)=hpr(k)*F1(k) + tau0*hjr(k)*F2(k)/sig^2;
    
end

hfpt=-hj0(1)/hjr(1); 

[dummy,P0,J0,r0]=EIF_r0_thy(DT,VT,Vth,Vre,tau0,tref,sig,E0,dV,Vl);

if w~=0
   hstr=hfpt/(1-hfpt); 
else
   hstr=pi*r0/dw0;          % this requires a dw to be correctly defined
end

hfsp=r0*(1+2*real(hstr));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,x]=ift(dw,w,hx);

nf=length(w);
dt=0.5; T1=-100; T2=500;
t=[T1:dt:T2]'; nt=length(t);
x=zeros(nt,1);

for s=1:nt
    x(s)=(1/(2*pi))*sum(dw.*exp(i*w*t(s)).*hx);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tisi,sisi,tstr,sstr,ffsp,sfsp]=sims(tau,sig,EL,Vth,Vre,VT,DT,tref);

nsp=20000; 
tsp=zeros(nsp,1); 

t=0; dt=0.1; 
noise=sqrt(2*dt/tau)*sig;
dtbytau=dt/tau;
tsp(1)=t; k=1;

refc=0;
Vold=EL;

while k<=nsp

    Vnew=Vold + (refc<=0)*((dtbytau)*(EL - Vold + DT*exp((Vold-VT)/DT)) + noise*randn); 
    
    if Vnew>Vth
        k=k+1; 
        tsp(k)=t;
        Vnew=Vre;
        refc=tref/dt;
    end

    t=t+dt;
    Vold=Vnew;
    refc=refc-1;
    
end

dt=2; % this is now the time step for the histograms of the spike data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the isi

isi=diff(tsp); nisi=length(isi);
%tisi=[0:dt:max(isi)]'; 
tisi=[0:dt:500]';
sisi=hist(isi,tisi)/(nisi*dt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the str

T=500;          % this is the range for time
cisi=cumsum(isi); cisimax=cisi(nisi);
s=max(find(cisi+T<cisimax));

str=[];

for k=1:s
    kstr=cisi(k+1:nisi)-cisi(k);
    pos=find(kstr<T);   
    str=[str;kstr(pos)];
end

nstr=length(str);
tstr=[0:dt:T]';
sstr=hist(str,tstr)/(s*dt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the fsp

fmaxHz=500;
T=400;
dt=1;

S=zeros(ceil(max(cisi/dt)),1);
cisi(1)=1;
S(ceil(cisi/dt),1)=1/dt;
s=ceil(T/dt); s1=1+s; s2=length(S)-s-1;

tAC=[-T:dt:T]'; nAC=length(tAC); 
AC=zeros(nAC,1);

for k=[-s:s]
    AC(k+s+1)=sum(S(s1+k:s2+k).*S(s1:s2))/(s2-s1+1);
end

% now fourier transform it...
nffsp=floor(2*T*fmaxHz/1000); 
ffsp=(1/(2*T))*[1:nffsp]';
wfsp=ffsp*2*pi;

sfsp=zeros(nffsp,1);

for k=1:nffsp
   sfsp(k)=sum(dt*exp(-i*wfsp(k)*tAC).*AC); 
end

sfsp=real(sfsp);   % cleans up and imaginary bits

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tt,VV]=sim_example(tau,sig,EL,Vth,Vre,VT,DT,tref,T,dt);

n=ceil(T/dt);
tt=dt*[1:n]';

VV=zeros(n,1); 
VV(1)=Vre;
refc=tref/dt;

noise=sqrt(2*dt/tau)*sig;
dtbytau=dt/tau;

for k=1:n-1

    VV(k+1)=VV(k) + (refc<=0)*((dt/tau)*(EL - VV(k) + DT*exp((VV(k)-VT)/DT)) + noise*randn); 
    
    if VV(k+1)>Vth
        VV(k)=Vth;
        VV(k+1)=Vre;
        refc=tref/dt;
    end
    
    refc=refc-1;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional spike dressing for the refractory period
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

L1=1/1; L2=1/2; 
ts=2;

A=(Vth-Vre)*L2/( L1*(1-exp(L2*ts)*exp(-L2*tref)) + L2*(1-exp(-L1*ts)));
B=(Vth-Vre)*L1/( L1*(1-exp(L2*ts)*exp(-L2*tref)) + L2*(1-exp(-L1*ts)));
a=Vth+A*exp(-L1*ts); 
b=Vre-B*exp(L2*ts)*exp(-L2*tref);

ss=ceil(ts/dt);
sref=ceil(tref/dt);
t=dt*[1:sref]';

Vsp1=a-A*exp(L1*(t-ts)); Vsp2=b+B*exp(L2*(ts-t));
Vsp=Vsp2; Vsp(1:ss)=Vsp1(1:ss);

pos=find((VV(1:n-1)==Vth)&(VV(2:n)==Vre));
ns=length(pos);

for k=1:ns
   VV(pos(k):pos(k)+sref-1)=Vsp; 
end
 
VV=VV(1:n);

end