function [V w St] = aEIF(t,I,C,gL,EL,Vt,Vp,Vr,Dt,tauw,a,b,initV,initw,pflag)

% [V w St] = aEIF(t,I,C,gL,EL,Vt,Vp,Vr,Dt,tauw,a,b,initV,initw,pflag)
%
% t     - time vector in sec. (deltaT might influence the results; smaller values earlier spikes; 0.0005)
% I     - injected current (in nA); 1. dim: time course; 2. dim: different neurons
% C     - membrane conductance (in nF)
% gL    - leak conductance (in nS)
% EL    - leak reversal (in mV), e.g., -70
% Vt    - threshold potential for spike (in mV)
% Vp    - peak potential at spike (in mV), e.g. 20
% Vr    - reset potential after spike (in mV), e.g. Vr = EL; Vr > Vt -> bursting
% Dt    - rise slope factor (in mV); sharpness of spike
% tauw  - adaptation time constant (~Ca-activated K current inactivation; in ms)
% a     - 'sub-threshold' adaptation conductance (in nS); large a = subthreshold oscillations
% b     - 'sra' current increment (in nA); large b = strong spike-frequency adaptation
% initV - initial value of membrane potential V (in mV; default = EL)
% initw - initial value of adaptation current w (in nA; default = 0)
% pflag - plotting 1 or 0
%
% V  - membrane potential (in mV)
% w  - adaptation current (in nA)
% St - spike times
%
% for example (tonic bursting):
% t = 0:0.0005:2; I = zeros([length(t) 1]); I(400:3601)=0.65;
% [V w] = aEIF(t,I,0.45,25,-70.6,-50.4,20,-47.4,2,144,4,0.0805,-70.6,0,1);
%
% for example (post-inhibitory rebound):
% t = 0:0.0005:2; I = zeros([length(t) 1]); I(200:1000)=-0.8;
% [V w] = aEIF(t,I,0.45,25,-60,-50.4,20,-60,2,720,80,0.0805,-70.6,0,1);
%
% for example (subthreshold oscillations and resonant firing):
% t = 0:0.0005:2; I = zeros([length(t) 1]); I([500:540 2200:2240 2400:2440])=1;
% [V w] = aEIF(t,I,0.576,6,-70.6,-50.4,20,-70.6,2,96,55,0.0805,-70.6,0,1);
%
% References:
% Brette R, Gerstner W (2005) Adaptive Exponential Integrate-and-Fire Model as an
%     Effective Description of Neuronal Activity. J Neurophysiol 94:3637-3642.
%
% https://wiki.brandeis.edu/twiki/bin/view/Labs/PMlab/AdExaEIF
% http://www.scholarpedia.org/article/Adaptive_exponential_integrate-and-fire_model
%
% Description: The scripts implements an adapting exponential integrate-and-fire 
% model (aEIF). Given an input current (I), the membrane potential (V) is calculated.
%
% ---------
%
%    This program 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.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% ------------------------------------------------------------------------------
% B. Herrmann, Email: herrmann.b@gmail.com, 2015-04-21

% convert to SI units
C     = C*1e-9;
gL    = gL*1e-9;
EL    = EL*1e-3;
Vt    = Vt*1e-3;
Vp    = Vp*1e-3;
Vr    = Vr*1e-3;
Dt    = Dt*1e-3;
a     = a*1e-9;
b     = b*1e-9;
tauw  = tauw*1e-3;
I     = I*1e-9;
initV = initV*1e-3;
initw = initw*1e-9;
dt    = t(2)-t(1);     % time step

% get number of neurons
nNeurons = size(I,2);

[V w Sb] = deal(zeros([length(t) nNeurons]));
V(1,:)   = initV;  % membrane potential vector
w(1,:)   = initw;  % adaptation current
for ii = 1 : length(t)-1
	% update membrane potential
	V(ii+1,:) = V(ii,:) + dt/C*(-gL*(V(ii,:)-EL) + gL*Dt*exp((V(ii,:)-Vt)/Dt) - w(ii,:) + I(ii,:));
	
	% check for spikes
	ixf = V(ii+1,:) > Vp;
	
	% update adaptation conductance
	w(ii+1,:) = w(ii,:) + dt/tauw*(a*(V(ii,:)-EL) - w(ii,:)) + ixf*b;
	
	% set spike and reset membrane potential
	V(ii,ixf)   = Vp;
	V(ii+1,ixf) = Vr;
	Sb(ii,ixf)  = 1; % spike bits
end

% get spike times
if nNeurons == 1
    St = find(Sb)*dt;
else
    St = cell([nNeurons 1]);
    for ii = 1: nNeurons
        St{ii} = find(Sb(:,ii))*dt;
    end
end

% convert back to original units
V = V*1e3;
w = w*1e9;
I = I*1e9;

% plot
if pflag
	figure, set(gcf,'Color',[1 1 1])
	subplot(6,1,1:3)
	plot(t,V,'-r')
	title('AdEx LIF neuron')
	ylabel('V (mV)')
	axis([t([1 end]) -100 30])
	subplot(6,1,4:5)
	plot(t,w,'-r')
	ylabel('w (nA)')
	axis([t([1 end]) min(w(:))-0.02 max(w(:))+0.02])
	subplot(6,1,6)
	plot(t,I,'-r')
	xlabel('t (s)'); ylabel('I (nA)')
	axis([t([1 end]) min(I(:))-0.02 max(I(:))+0.02])
end


