%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Rate response for a LIF neuron receiving Poissonian white noise.
%
% This script runs in MATLAB (and hopefully OCTAVE) and will generate 
% panels A,B,C and D of figure 2 in the paper:
%
% Firing-rate response of a neuron receiving excitatory and inhibitory synaptic shot noise
% Richardson MJE and Swarbrick R
% Physical Review Letters 105:article-no 178102 (2010)
%
% Please reference this publication in any document that uses the code.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COPYRIGHT: Magnus Richardson (2011).
% 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 RandS2010fig2ABCD

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

vth=10; vre=5; tau=20;

f1=1/K; f2=1000000/K; nf=49;
fkHz=exp(linspace(log(f1),log(f2),nf))';
fHz=K*fkHz; w=2*pi*fkHz;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeA=2; aiA=-0.1; sig0A=4; mu0A=5;
Re0A=(sig0A^2-mu0A*aiA)/(tau*aeA*(aeA-aiA)); Re1A=0.02;
Ri0A=(sig0A^2-mu0A*aeA)/(tau*aiA*(aiA-aeA)); Ri1A=0.6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeB=0.1; aiB=-0.1; sig0B=4; mu0B=5;
Re0B=(sig0B^2-mu0B*aiB)/(tau*aeB*(aeB-aiB)); Re1B=0.3; 
Ri0B=(sig0B^2-mu0B*aeB)/(tau*aiB*(aiB-aeB)); Ri1B=0.3; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeC=0.1; aiC=-2; sig0C=4; mu0C=5;
Re0C=(sig0C^2-mu0C*aiC)/(tau*aeC*(aeC-aiC)); Re1C=0.3;
Ri0C=(sig0C^2-mu0C*aeC)/(tau*aiC*(aiC-aeC)); Ri1C=0.031;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D example (diffusion case)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mue1=tau*Re1B*aeB;
mui1=tau*Ri1B*aiB;
sige1=sqrt(tau*Re1B*aeB^2);
sigi1=sqrt(tau*Ri1B*aiB^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[r0A,are1A,thre1A,ari1A,thri1A,re1A,ri1A]=re1ri1_exact1(tau,vre,vth,Re0A,Ri0A,aeA,aiA,Re1A,Ri1A,w);
[r0B,are1B,thre1B,ari1B,thri1B,re1B,ri1B]=re1ri1_exact1(tau,vre,vth,Re0B,Ri0B,aeB,aiB,Re1B,Ri1B,w);
[r0C,are1C,thre1C,ari1C,thri1C,re1C,ri1C]=re1ri1_exact1(tau,vre,vth,Re0C,Ri0C,aeC,aiC,Re1C,Ri1C,w);
[r0D,are1D,thre1D,ari1D,thri1D,re1D,ri1D]=re1ri1_gauss(tau,vre,vth,sig0B,mu0B,mue1,mui1,sige1,sigi1,w); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,1); 
loglog(fHz,K*are1A/Re1A,'g',fHz,K*are1B/Re1B,'b',fHz,K*are1C/Re1C,'r',fHz,K*are1D/Re1B,'k'); 
axis([min(fHz) max(fHz) 0.1 1000]);
xlabel('Freq [Hz]'); ylabel('Rate/Re1'); title('Excitation');

subplot(2,2,3); 
semilogx(fHz,r2d*thre1A,'g',fHz,r2d*thre1B,'b',fHz,r2d*thre1C,'r',fHz,r2d*thre1D,'k'); 
axis([min(fHz) max(fHz) -60 15]);
xlabel('Freq [Hz]'); ylabel('Phase (deg)');

subplot(2,2,2); 
loglog(fHz,K*ari1A/Ri1A,'g',fHz,K*ari1B/Ri1B,'b',fHz,K*ari1C/Ri1C,'r',fHz,K*ari1D/Ri1B,'k'); 
axis([min(fHz) max(fHz) 0.001 1000]);
xlabel('Freq [Hz]'); ylabel('Rate/Ri1'); title('Inhibition');

subplot(2,2,4); 
semilogx(fHz,r2d*thri1A,'g',fHz,r2d*thri1B,'b',fHz,r2d*thri1C,'r',fHz,r2d*thri1D,'k'); 
axis([min(fHz) max(fHz) -15 180]);
xlabel('Freq [Hz]'); ylabel('Phase (deg)');

drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% excitatory rate modulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,are1,thre1,ari1,thri1,re1,ri1]=re1ri1_exact1(tau,vre,vth,Re0,Ri0,ae,ai,Re1,Ri1,w)

K=1000;

% get the steady state
ds=(1/ae)/20000; s=[ds/2:ds:1/ae-ds/2];
A=exp(Re0*tau*log(1-ae*s)+Ri0*tau*log(1-ai*s));
B=(exp(vth*s)./(1-ae*s)-exp(vre*s))./s;
m=find(isnan(A.*B)); A(m)=0; B(m)=0;
I=sum(A.*B*ds); r0=1./(tau*I);

% now the dynamics
be=Re0*tau; bi=Ri0*tau; a=ai/ae;
nw=length(w); r1=zeros(nw,1);
fHz=w*K/(2*pi);

x1=0.0000001; x2=200; lx=500; 
x=exp(linspace(log(x1),log(x2),lx))';
logx=log(x); dlogx=logx(2)-logx(1); logx=[logx;logx(end)+dlogx]-dlogx/2;
dx=diff(exp(logx)); % makes the dx well spaced in logspace
Iy=zeros(lx,1);

re1=zeros(nw,1);
ri1=zeros(nw,1);

m1=0.0000001; m2=200; lm=500; 
m=exp(linspace(log(m1),log(m2),lm))';
logm=log(m); dlogm=logm(2)-logm(1); logm=[logm;logm(end)+dlogm]-dlogm/2;
dm=diff(exp(logm)); % makes the dx well spaced in logspace

for n=1:nw

L=w(n)*tau;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the saddle point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C=[(1+1i*L),-((1+a)*(1+i*L)+be+a*bi),a*(1+i*L+be+bi)];
R=roots(C); zstars=log(R); s=find(imag(zstars)<0&real(zstars)>0);

zstar=zstars(s); rstar=abs(zstar); angstar=angle(zstar);
eiangstar=exp(1i*angstar); degstar=angstar*360/(2*pi);

z=x*eiangstar; 
dz=dx*eiangstar; 

Dpsi=(z-zstar)*(1+1i*L)-be*log((1-exp(-z))/(1-exp(-zstar)))-bi*log((1-a*exp(-z))/(1-a*exp(-zstar)));
AbyAstar=exp(-Dpsi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% excitation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
for k=1:lx
    Iy(k,1)=(1/(1i*L+1))*sum(dm.*exp(-m)./(1-exp(-m/(1i*L+1)-z(k))));
end

% the numerator
s=exp(-z)/ae; 
f=(exp(s*vth)./(1-ae*s)-exp(s*vre));
fArIy=f.*AbyAstar.*Iy; 
F=sum(dz.*fArIy);

% now the denominator
s=exp(-z)/ae; 
g=((exp(s*vth)./(1-ae*s)-exp(s*vre)))./(ae*s);
gAr=g.*AbyAstar; 
G=sum(dz.*gAr);

re1(n)=Re1*tau*r0*F/G;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inhibition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:lx 
   Iy(k,1)=(1/(1+1i*L))*sum(dm.*exp(-m)./(1-a*exp(-m/(1+1i*L)-z(k))));
end

% the numerator
s=exp(-z)/ae; 
f=(exp(s*vth)./(1-ae*s)-exp(s*vre));
fArIy=f.*AbyAstar.*Iy; 
F=sum(dz.*fArIy);

% now the denominator
s=exp(-z)/ae; 
g=((exp(s*vth)./(1-ae*s)-exp(s*vre)))./(ae*s);
gAr=g.*AbyAstar; 
G=sum(dz.*gAr);

ri1(n)=Ri1*tau*r0*(ai/ae)*F/G;

end

are1=abs(re1); thre1=angle(re1);
ari1=abs(ri1); thri1=angle(ri1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% steady-state rate from the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,are1,thre1,ari1,thri1,re1,ri1]=re1ri1_gauss(tau,vre,vth,sig0,mu0,mue1,mui1,sige1,sigi1,w)

% calculate the steady state
yh=200; dy=0.0001; 
y=[dy/2:dy:yh-dy/2]';
yth=(vth-mu0)/sig0;
yre=(vre-mu0)/sig0;
A=exp(-y.^2/2).*(exp(y*yth)-exp(y*yre))./y;
r0=1/(tau*dy*sum(A));

nw=length(w);

da=0.0001; a=[da/2:da:1-da/2]';
db=0.0001; b=[db/2:db:100-db/2]';

psia=1i*log(1i*a/(1+1i))+(a.^2+2i)/4;
psib=1i*log((b+1i)/(1+1i))-((b+1i).^2-2i)/4;

re1=zeros(nw,1);
ri1=zeros(nw,1);

for k=1:nw
    
    wt=w(k)*tau;
    
    % for a int (f numerator for rho, g the denominator)
    Za=1i*exp(wt*psia); 
    y=1i*a*sqrt(wt/2);
    f0a=(exp(y*yth)-exp(y*yre))./y; 
    f1a=(exp(y*yth)-exp(y*yre)); 
    f2a=(exp(y*yth)-exp(y*yre)).*y; 
    F0a=Za.*f0a; F0a(isnan(F0a))=0;
    F1a=Za.*f1a; F1a(isnan(F1a))=0;
    F2a=Za.*f2a; F2a(isnan(F2a))=0;
     
    % for b int (f numerator for rho, g denominator)
    Zb=exp(wt*psib); 
    y=(b+1i)*sqrt(wt/2); 
    f0b=(exp(y*yth)-exp(y*yre))./y; 
    f1b=(exp(y*yth)-exp(y*yre)); 
    f2b=(exp(y*yth)-exp(y*yre)).*y; 
    F0b=Zb.*f0b; F0b(isnan(F0b))=0;
    F1b=Zb.*f1b; F1b(isnan(F1b))=0;
    F2b=Zb.*f2b; F2b(isnan(F2b))=0;
    
    F0=sum(F0a*da)+sum(F0b*db); 
    F1=sum(F1a*da)+sum(F1b*db); 
    F2=sum(F2a*da)+sum(F2b*db); 
    
    I1=F1/F0;
    I2=F2/F0;
       
    re1(k)=r0*(I1*mue1/(sig0*(1+1i*wt)) + I2*sige1^2/(sig0^2*(2+1i*wt)));
    ri1(k)=r0*(I1*mui1/(sig0*(1+1i*wt)) + I2*sigi1^2/(sig0^2*(2+1i*wt)));
    
end

are1=abs(re1); thre1=angle(re1);
ari1=abs(ri1); thri1=angle(ri1);