%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This file contains 4 independent functions that should run in either the
% OCTAVE or MATLAB environments. See text next to each function for a
% more complete description.
%
% LIF0
% generates the steady-state rate and density for the Leaky IF model
%
% LIF1
% generates the modulated rate for the Leaky IF model
%
% EIF0
% generates the steady-state rate and density for the Exponential IF model
%
% EIF1
% generates the modulated rate for the Exponential IF model
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The following publications should be referenced in any
% document that makes use of the Threshold Integration code:
%
% Firing-rate response of linear and nonlinear integrate-and-fire neurons
% to modulated current-based and conductance-based synaptic drive
% Richardson MJE, Physical Review E 76 article-no 021919 (2007)
%
% Spike-train spectra and network response functions for
% non-linear integrate-and-fire neurons
% Richardson MJE, Biological Cybernetics 99: 381-392 (2008).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COPYRIGHT: Magnus Richardson (2009).
% 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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LIF0
% Steady-state rate and density for the Leaky Integrate-and-Fire model
% using Threshold Integration.
% Runs in either the OCTAVE or MATLAB environment.
%
% The voltage obeys the equation:
% tau*dVdt = E0 - V + sig*sqrt(2*tau)*xi(t)
%
% where xi(t) is Gaussian white noise with =delta(t-t')
% such that integral(t->t+dt) xi(t)dt =sqrt(dt)*randn
%
% The threshold is at Vth with a reset at Vre.
%
% input: expected units are: tau [ms], sig,E0,Vth and Vre all in [mV]
% (example values might be tau=20; sig=5; E0=-50; Vth=-50; Vre=-60;)
%
% output: vectors V [mV] range of voltage
% P0 [1/mV] probability density
% J0 [kHz] probability flux (a piece-wise constant)
% scalar r0 [kHz] steady state firing rate
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,P0,J0,r0]=LIF0(tau,sig,E0,Vth,Vre)
% set up the lattice for the discretized integration
dV=0.01;
Vl=-100; % the lower bound for the voltage range
V=[Vl:dV:Vth]';
n=length(V);
kre=find(V==Vre); % NB Vre must fall on a lattice point!
% Will use the modified Euler method (Phys. Rev. E 2007, Eqs A1-A6)
G=(V-E0)/sig^2; % the only part specific to the Leaky IF
A=exp(dV*G);
B=(A-1)./(sig^2*dV*G);
B(find(G==0))=1/sig^2;
% set up the vectors for the scaled current and probability density
j0=zeros(n,1);
p0=zeros(n,1);
j0(n)=1; % initial conditions at V=Vth. NB p0(n)=0 already.
for k=n:-1:2 % integrate backwards from threshold
j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*A(k) + dV*B(k)*tau*j0(k);
end
r0=1/(dV*sum(p0)); % steady-state firing rate (in kHz)
P0=r0*p0; % correctly normalised density and current
J0=r0*j0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LIF1
% Leaky Integrate-and-Fire model with modulated parameters
% using Threshold Integration.
% Runs in either the OCTAVE or MATLAB environment.
%
% The code is written for current modulation (E) but can be easily modified
% for conductance or noise modulation.
%
% The voltage obeys the equation:
% tau*dVdt = E0 + E1*cos(w*t) - V + sig*sqrt(2*tau)*xi(t)
%
% where xi(t) is Gaussian white noise with =delta(t-t')
% such that integral(t->t+dt) xi(t)dt =sqrt(dt)*randn
%
% The threshold is at Vth with a reset at Vre.
%
% The firing rate response is calculated to first order
% r(t)=r0 + abs(r1)*cos(w*t+phase(r1))
%
% input: expected units are: tau [ms], sig,E0,E1,Vth and Vre all in [mV]
% and fkHz is the modulated frequency f=w/2*pi in kHz
%
% output: r0 [kHz] steady state firing rate
% r1 [kHz] rate response (complex number)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,r1]=LIF1(tau,sig,E0,Vth,Vre,E1,fkHz)
% set up the voltage range
Vl=-100;
dV=0.01;
V=[Vl:dV:Vth]';
n=length(V);
kre=find(V==Vre);
% Will use the modified Euler method (Phys. Rev. E 2007, Eqs A1-A6)
G=(V-E0)/sig^2; % the only part specific to the Leaky IF
A=exp(dV*G);
B=(A-1)./(sig^2*dV*G);
B(find(G==0))=1/sig^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first need the steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j0=zeros(n,1); j0(n)=1;
p0=zeros(n,1);
for k=n:-1:2 % integrate backwards from threshold
j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*A(k) + dV*B(k)*tau*j0(k);
end
r0=1/(dV*sum(p0)); % steady-state firing rate (in kHz)
P0=r0*p0; % normalised density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the response to current (E) modulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=2*pi*fkHz;
jh1=zeros(n,1); jh1(n)=1;
ph1=zeros(n,1);
jhE=zeros(n,1);
phE=zeros(n,1);
for k=n:-1:2
jh1(k-1)=jh1(k) + dV*i*w*ph1(k) - (k==kre+1);
ph1(k-1)=ph1(k)*A(k) + dV*B(k)*tau*jh1(k);
% the 2nd equation in the following pair is specific to E modulation
% if conductance or noise modulation is required then
% it is only the last (inhomogeneous) term containing P0 that must be
% altered - please see Eq 33 of Phys. Rev. E (2007) paper for details
jhE(k-1)=jhE(k) + dV*i*w*phE(k);
phE(k-1)=phE(k)*A(k) + dV*B(k)*(tau*jhE(k) - E1*P0(k));
end
r1=-jhE(1)/jh1(1); % because Jh(Vl)=0 => jhe(1) + r1*jh1(1)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EIF0
% Steady-state rate and density for the
% Exponential Integrate-and-Fire model using Threshold Integration.
% Runs in either the OCTAVE or MATLAB environment.
%
% The voltage obeys the equation:
% tau*dVdt = E0 - V + DT*exp((V-VT)/DT) + sig*sqrt(2*tau)*xi(t)
%
% where xi(t) is Gaussian white noise with =delta(t-t')
% such that integral(t->t+dt) xi(t)dt =sqrt(dt)*randn
%
% The threshold is at Vth with a reset at Vre.
%
% input: expected units are: tau [ms], sig,E0,Vth,Vre,VT,DT all in [mV]
% example values might be tau=20; sig=5; E0=-50;
% Vth=0; Vre=-60; DT=3; VT=-53;
%
% output: vectors V [mV] range of voltage
% P0 [1/mV] probability density
% J0 [kHz] probability flux (a piece-wise constant)
% scalar r0 [kHz] steady state firing rate
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,P0,J0,r0]=EIF0(tau,sig,E0,Vth,Vre,VT,DT)
% set up the lattice for the discretized integration
dV=0.01;
Vl=-100; % the lower bound for the voltage range
V=[Vl:dV:Vth]';
n=length(V);
kre=find(V==Vre); % NB Vre must fall on a lattice point!
% Will use the modified Euler method (Phys. Rev. E 2007, Eqs A1-A6)
G=(V-E0-DT*exp((V-VT)/DT))/sig^2; % only bit specific to the EIF
A=exp(dV*G);
B=(A-1)./(sig^2*dV*G);
B(find(G==0))=1/sig^2;
% set up the vectors for the scaled current and probability density
j0=zeros(n,1);
p0=zeros(n,1);
j0(n)=1; % initial conditions at V=Vth. NB p0(n)=0 already.
for k=n:-1:2 % integrate backwards from threshold
j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*A(k) + dV*B(k)*tau*j0(k);
end
r0=1/(dV*sum(p0)); % steady-state firing rate (in kHz)
P0=r0*p0; % correctly normalised density and current
J0=r0*j0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EIF1
% Exponential Integrate-and-Fire model with modulated parameters
% using Threshold Integration.
% Runs in either the OCTAVE or MATLAB environment.
%
% The code is written for current modulation (E) but can be easily modified
% for conductance or noise modulation.
%
% The voltage obeys the equation:
% tau*dVdt=E0 + E1*cos(w*t) - V + DT*exp((V-VT)/DT) + sig*sqrt(2*tau)*xi(t)
%
% where xi(t) is Gaussian white noise with =delta(t-t')
% such that integral(t->t+dt) xi(t)dt =sqrt(dt)*randn
%
% The threshold is at Vth with a reset at Vre.
%
% The firing rate response is calculated to first order
% r(t)=r0 + abs(r1)*cos(w*t+phase(r1))
%
% input: expected units are: tau [ms], sig,E0,E1,Vth,Vre,DT,VT all in [mV]
% and fkHz is the modulated frequency f=w/2*pi in kHz
%
% output: r0 [kHz] steady state firing rate
% r1 [kHz] rate response (complex number)
%
% NOTE: The following publications should be referenced in any
% document that makes use of the code:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,r1]=EIF1(tau,sig,E0,Vth,Vre,VT,DT,E1,fkHz)
% set up the voltage range
Vl=-100;
dV=0.01;
V=[Vl:dV:Vth]';
n=length(V);
kre=find(V==Vre);
% Will use the modified Euler method (Phys. Rev. E 2007, Eqs A1-A6)
G=(V-E0-DT*exp((V-VT)/DT))/sig^2; % only bit specific to the EIF
A=exp(dV*G);
B=(A-1)./(sig^2*dV*G);
B(find(G==0))=1/sig^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first need the steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j0=zeros(n,1); j0(n)=1;
p0=zeros(n,1);
for k=n:-1:2 % integrate backwards from threshold
j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*A(k) + dV*B(k)*tau*j0(k);
end
r0=1/(dV*sum(p0)); % steady-state firing rate (in kHz)
P0=r0*p0; % normalised density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the response to current (E) modulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=2*pi*fkHz;
jh1=zeros(n,1); jh1(n)=1;
ph1=zeros(n,1);
jhE=zeros(n,1);
phE=zeros(n,1);
for k=n:-1:2
jh1(k-1)=jh1(k) + dV*i*w*ph1(k) - (k==kre+1);
ph1(k-1)=ph1(k)*A(k) + dV*B(k)*tau*jh1(k);
% the 2nd equation in the following pair is specific to E modulation
% if conductance or noise modulation is required then
% it is only the last (inhomogeneous) term containing P0 that must be
% altered - please see Eq 33 of Phys. Rev. E (2007) paper for details
jhE(k-1)=jhE(k) + dV*i*w*phE(k);
phE(k-1)=phE(k)*A(k) + dV*B(k)*(tau*jhE(k) - E1*P0(k));
end
r1=-jhE(1)/jh1(1); % because Jh(Vl)=0 => jhe(1) + r1*jh1(1)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%