%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Steady-state rate for a LIF neuron receiving Poissonian white noise.
%
% This script runs in MATLAB (and hopefully OCTAVE) and will generate
% panel A of figure 1 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 RandS2010fig1A
clf
K=1000;
tau=20; vre=5; vth=10; vlb=-20;
mu0lb=-10; mu0ub=10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Case A. Large excitatory jumps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeA=2; aiA=-0.1; sig0A=4;
mu0A=[max(mu0lb,sig0A^2/aiA):0.5:min(sig0A^2/aeA,mu0ub)]'; nA=length(mu0A);
Re0A=(sig0A^2-mu0A*aiA)/(tau*aeA*(aeA-aiA));
Ri0A=(sig0A^2-mu0A*aeA)/(tau*aiA*(aiA-aeA));
r0thA=zeros(nA,1);
for k=1:nA
r0thA(k,1)=r0exact(tau,vre,vth,Re0A(k),Ri0A(k),aeA,aiA);
end
'case A run complete'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% case B. Diffusion case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeB=0.1; aiB=-0.1; sig0B=4;
mu0B=[max(mu0lb,sig0B^2/aiB):0.5:min(sig0B^2/aeB,mu0ub)]'; nB=length(mu0B);
Re0B=(sig0B^2-mu0B*aiB)/(tau*aeB*(aeB-aiB));
Ri0B=(sig0B^2-mu0B*aeB)/(tau*aiB*(aiB-aeB));
r0thB=zeros(nB,1);
r0thD=zeros(nB,1); % for the diffusion approximation
for k=1:nB
r0thB(k,1)=r0exact(tau,vre,vth,Re0B(k),Ri0B(k),aeB,aiB);
r0thD(k,1)=r0exactdiff(tau,vre,vth,Re0B(k),Ri0B(k),aeB,aiB);
end
'case B run complete'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% case C. Large inhibitory jumps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
aeC=0.1; aiC=-2; sig0C=4;
mu0C=[max(mu0lb,sig0C^2/aiC):0.5:min(sig0C^2/aeC,mu0ub)]'; nC=length(mu0C);
Re0C=(sig0C^2-mu0C*aiC)/(tau*aeC*(aeC-aiC));
Ri0C=(sig0C^2-mu0C*aeC)/(tau*aiC*(aiC-aeC));
r0thC=zeros(nC,1);
for k=1:nC
r0thC(k,1)=r0exact(tau,vre,vth,Re0C(k),Ri0C(k),aeC,aiC);
end
'case C run complete'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the figures for theory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(mu0A,K*r0thA,'g',mu0B,K*r0thB,'b',mu0B,K*r0thD,'k:',mu0C,K*r0thC,'r');
ylabel('Rate r0 [Hz]');
xlabel('mu_0 [mV]');
axis([-10 10 0 25]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% steady-state rate from the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r0=r0exact(tau,vre,vth,Re0,Ri0,ae,ai)
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% steady-state rate from the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r0=r0exactdiff(tau,vre,vth,Re0,Ri0,ae,ai)
mu0=ae*Re0*tau + ai*Ri0*tau;
sig0=sqrt(ae^2*Re0*tau + ai^2*Ri0*tau);
yh=100; dy=0.01; 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;
IA=dy*sum(A);
r0=1/(tau*IA);