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