%2017-09-18

% Integrate & Fire model with spike-rate adaptation
% Based on the passive membrane equation
% Numerical solution using the modified Euler method (Runge-Kutta, order 2)

% Ref: Honeycutt R, "Stochastic Runge-Kutta algorithm I. White noise",
% Phys. Rev. A, 45:600-603



clearvars;
close all;

% [C] = muF/cm^2
% [V] = mV
% [G] = mS/cm^2
% [Iapp] = muA/cm^2

% Biophysical parameters

C = 1;
El = -60;
Ek = -80;
Gl = 0.1;
Gk = 0.01;
Iapp = 1.01;
D = 0;

% Threshold and reset values for V

Vth = -50;
Vrst = -65;

% spike-rate adaptation parameters

tausra = 100;
Dtawsra = 0.01;

% Time definitions

Tmax = 4000;
dt = 0.1;
t = 0:dt:Tmax;

% Square wave function

ti = dt;
tf = 10000;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;

% Initial conditions

V = zeros(1,length(t));
w = zeros(1,length(t));
V(1) = El;
w(1) = 0;

% Computation of the solution

for j=1:length(t)-1
    kv1 = (-Gl*(V(j)-El)-Gk*w(j)*(V(j)-Ek)+Iapp*H(j))/C;
    kw1 = -w(j)/tausra;
    av = V(j)+kv1*dt;
    aw = w(j)+kw1*dt;
    kv2 = (-Gl*(av-El)-Gk*aw*(av-Ek)+Iapp*H(j+1))/C;
    kw2 = -aw/tausra;
    V(j+1) = V(j) + (kv1+kv2)*dt/2;
    w(j+1) = w(j)+(kw1+kw2)*dt/2;
    if V(j+1)>=Vth
       V(j)=60;
       V(j+1)=Vrst;
       w(j+1)=w(j+1)+Dtawsra;
   end
end

% Initial conditions (white noise)

tspk = zeros(1);
ISI = zeros(1);
Vn = zeros(1,length(t));
wn = zeros(1,length(t));
Vn(1) = V(1);
wn(1) = 0;

% Computation of the solution (white noise)

spkcnt = 0;
for j=1:length(t)-1
    eta = randn;
    Vnaux = Vn(j)+sqrt(2*D*dt)*eta;
    kv1 = (-Gl*(Vnaux-El)-Gk*wn(j)*(Vn(j)-Ek)+Iapp*H(j))/C;
    kw1 = -wn(j)/tausra;
    av = Vn(j)+kv1*dt;
    av = av + sqrt(2*D*dt)*eta;
    aw = wn(j)+kw1*dt;
    kv2 = (-Gl*(av-El)-Gk*aw*(av-Ek)+Iapp*H(j+1))/C;
    kw2 = -aw/tausra;
    Vn(j+1) = Vn(j) + (kv1+kv2)*dt/2;
    Vn(j+1) = Vn(j+1) + sqrt(2*D*dt)*eta;
    wn(j+1) = wn(j)+(kw1+kw2)*dt/2;
    if Vn(j+1)>=Vth
       Vn(j)=60;
       Vn(j+1)=Vrst;
       spkcnt = spkcnt+1;
       tspk(spkcnt)=t(j);
       wn(j+1)=wn(j+1)+Dtawsra;
   end
end

% Computation of the inter-spike intervals (ISI)

for k=1:length(tspk)-1;
    ISI(k) = tspk(k+1)-tspk(k);
end

% Graphs

figure
hold on
plot(t,V,'b','linewidth',2);
if D > 0
    plot(t,Vn,'r','linewidth',2);
end
axis([0 Tmax -80 80]);
set(gca,'fontsize',20);
xlabel('t')
ylabel('V');

if D>0
    figure
    hist(ISI,20);
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end

