%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 <http://www.gnu.org/licenses/> 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 <xi(t)xi(t')>=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 <xi(t)xi(t')>=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 <xi(t)xi(t')>=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 <xi(t)xi(t')>=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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%