%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 <http://www.gnu.org/licenses/> 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);