%2020-01-03

% Quadratic 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;
Gqif = 0.1;
Vr = -65;
Vc = -55;
Iapp = 2.6;
D = 0;

% Threshold and reset values for V

Vth = -50;
Vrst = -70;

% 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 = (Gqif*(V(j)-Vr)*(V(j)-Vc)+Iapp)/C;
    av = V(j)+kv1*dt;
    kv2 = (Gqif*(av-Vr)*(av-Vc)+Iapp)/C;
    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');

% Phase-space diagram

vv=-100:0.001:100;
F = Gqif*(vv-Vr).*(vv-Vc)+Iapp;

if Gqif>0
    Vfp = roots([1 -(Vr+Vc) Vr*Vc+Iapp/Gqif]);
end
       
        
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 Gqif>0 && imag(Vfp(1)) == 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 = (Gqif*(Vn(j)-Vr)*(Vn(j)-Vc)+Iapp)/C;
        av = Vn(j)+kv1*dt;
        av = av + sqrt(2*D*dt)*eta;
        kv2 = (Gqif*(av-Vr)*(av-Vc)+Iapp)/C;
        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
