%2020-01-03

% Exponential Integrate & Fire model
% 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

% Parameters

C = 1;
Geif = 0.1;
Vr = -65;
Vrh = -55;
DeltaT = 2;
Iapp = 2;
D = 0;

% Threshold and reset values for V

Vth = -30;
Vrst = -60;

% Time definitions

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

% Initial conditions

V = zeros(1,length(t));
V(1) = Vrst;

% Computation of the solution

tspk = zeros(1);

spkcnt = 0;
for j=1:length(t)-1
    kv1 = -Geif*(V(j)-Vr)+Geif*DeltaT*exp((V(j)-Vrh)/DeltaT)+Iapp;
    av = V(j)+kv1*dt;
    kv2 = -Geif*(av-Vr)+Geif*DeltaT*exp((av-Vrh)/DeltaT)+Iapp;
    V(j+1) = V(j) + (kv1+kv2)*dt/2;
    if V(j+1)>=Vth
       V(j)=60;
       V(j+1)=Vrst;
       spkcnt = spkcnt+1;
       tspk(spkcnt)=t(j);
   end
end

% Computation of the inter-spike intervals (ISI)

ISI = diff(tspk);

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

% Computation of the fixed-point

vv=-100:0.001:100;
Feif =@(V) -Geif*(V-Vr)+Geif*DeltaT*exp((V-Vrh)/DeltaT)+Iapp;
Vfp = zeros(1);

fpcnt = 0;
for j=1:length(vv)
    if Feif(vv(j))<0 && fpcnt == 0
        fpcnt = 1;
        Vfp(fpcnt) = (vv(j)+vv(j-1))/2;        
    elseif Feif(vv(j))>0 && fpcnt == 1
        fpcnt = 2;
        Vfp(fpcnt) = (vv(j)+vv(j-1))/2;        
    end
end


% Phase-space diagram


F = Feif(vv);

figure
hold on
plot(vv,F,'r','linewidth',2);
plot([-100 100],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
axis([-80 -20 -10 20]);
set(gca,'fontsize',20);
xlabel('V  [mV]');
ylabel('dV / dt');
plot([Vth Vth],[-20 20],'--','Color',[.7 .7 .7],'linewidth',2);
plot([Vrst Vrst],[-20 20],'--','Color',[.7 .7 .7],'linewidth',2);
if fpcnt > 0
    plot(Vfp(1),0,'ok','linewidth',2);
    plot(Vfp(1),0,'xk','linewidth',2);
    plot(Vfp(2),0,'ok','linewidth',2); 
end
text(Vth,-9,'V_{th}','fontsize',20)
text(Vrst,-9,'V_{reset}','fontsize',20)


if D > 0 

    Tmax = 500000;
    dt = 0.1;
    t = 0:dt:Tmax;
    
    % Initial conditions (white noise)

    tspkn = zeros(1);    
    Vn = zeros(1,length(t));
    Vn(1) = V(1);

    % Computation of the solution (white noise)

    spkcntn = 0;
    for j=1:length(t)-1
        eta = randn;
        kv1 = -Geif*(Vn(j)-Vr)+Geif*DeltaT*exp((Vn(j)-Vrh)/DeltaT)+Iapp;
        av = Vn(j)+kv1*dt;
        av = av + sqrt(2*D*dt)*eta;
        kv2 = -Geif*(av-Vr)+Geif*DeltaT*exp((av-Vrh)/DeltaT)+Iapp;
        Vn(j+1) = Vn(j) + (kv1+kv2)*dt/2;
        Vn(j+1) = Vn(j+1) + sqrt(2*D*dt)*eta;
        if Vn(j+1)>=Vth
           Vn(j)=60;
           Vn(j+1)=Vrst;
           spkcntn = spkcntn+1;
           tspkn(spkcntn)=t(j);
       end
    end
    
    % Computation of the inter-spike intervals (ISI)

    ISIn = diff(tspkn);
   

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

    Nhist = 100;
    figure
    hold on
    hist(ISIn,Nhist); 
    if tspk(end)>0
        plot([ISI(end) ISI(end)],[0 max(hist(ISIn,Nhist))],'-r','linewidth',2);
    end
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
   
end


