% Leaky integrate-and-fire Neuron (LIF), solved with Forward Euler (FE),
% Backward Euler (BE) and Crank-Nicolson (CN). A constant current is applied,
% such that an analytical solution is known. This allows to measure convergence
% rates of the different schemes (FE, BE: 1; CN:2).
%
% === parameter ===
%
% Parameter from Dayan-Abbott: Theoretical Neuroscience
% (But the current values are a rough guess estimated from an equilibrium
% condition [which here is I = 1.5 nA with constant current]
% to get meaningful values)
%
% === analytical solution ===
%
% In this script a constant current I is assumed to stimulate the Neuron. An
% analytical solution is known fo that case:
%
% vana = vsteady + [vinit - vsteady] * exp(-t / taum)
%
% in the interspiking intervals where vsteady = vleak + I*R. The cell is
% initialized (and reset after spike) to its resting potential vrest = vleak,
% thus we have vinit = vleak here, and t is relative to the last reset time.

%------------------------------------------------------------------
% create a structure for the leaky integrate and fire neuron (LIF)
%------------------------------------------------------------------
lif.R = 10.;                 % resistancy [MOhm]
lif.taum = 10.;              % time constant [ms]
lif.vreset = -65.;           % reset value for potential (resting pot) [mV]
lif.vleak = -65.;            % leak potential [mV] (= resting pot.)
lif.vpeak = -50.;            % spike detection threshold [mV]
lif.tdet = 0.;               % time of spike detection [ms]
lif.vdet = 0.;               % potential at tdet [mV]
lif.T = 10.;                 % simulation time [ms]
lif.dt = 0.1;                % time step size [ms]
lif.n = round(lif.T/lif.dt); % number of simulation steps
lif.I = 2.0;                 % constant current in nA

%---------------------
% analytical solution
%---------------------
function ret = anaSol(lif, t)
   vsteady = lif.vleak + lif.R * lif.I;
   ret = vsteady + (lif.vleak - vsteady) * exp(-t / lif.taum);
end;

%---------------------------------
% solve with Forward Euler method
%---------------------------------
function ret = solveFE(lif)
    t = [];
    tlocal = 0.;
    v = lif.vleak;
    for i = 1 : lif.n-1
        % FE
        t(i)  = i * lif.dt;
        v = v + lif.dt/lif.taum * (lif.R * lif.I - (v - lif.vleak));

        % spike detected !
        if v >= lif.vpeak
            v  = lif.vreset; % membrane potential reset
            tlocal = t(i);
        end;
        
        % evaluate analytical solution
        vana = anaSol(lif, -tlocal + t(i));
        errV = v - vana;

        % save data in vectors for plotting
        ii(i)  = lif.I;
        vv(i)  = v;
        err(i) = errV;

        % print to standard cout
        %printf("%g\n", v);
    end;

    errV

    %% plot  solution
    %figure;
    %plot(t, vv, t, ii); grid;
    %xlabel('t [ms]');
    %ylabel('V [mV] and I [nA]');
    %legend('Potential', 'Applied current');
    %replot
    
    %% export
    %print ("iaf_leaky-potential.pdf", "-color", "-dpdf");
end;

%----------------------------------
% solve with Backward Euler method
%----------------------------------
function ret = solveBE(lif)
    t = [];
    tlocal = 0.;
    v = lif.vleak;
    for i = 1 : lif.n-1
        % BE
        t(i)  = i * lif.dt;
        v = 1.0 / (1.0 + lif.dt / lif.taum) * (v + lif.dt/lif.taum * (lif.R * lif.I + lif.vleak));

        % spike detected !
        if v >= lif.vpeak
            v  = lif.vreset; % membrane potential reset
            tlocal = t(i)
        end;
        
        % evaluate analytical solution
        vana = anaSol(lif, -tlocal + t(i));
        errV = v - vana;

        % save data in vectors for plotting
        ii(i)  = lif.I;
        vv(i)  = v;
        err(i) = errV;

        % print to standard cout
        %printf("%f\n", v);
    end;

    errV

    %% plot  solution
    %figure;
    %plot(t, vv, t, ii); grid;
    %xlabel('t [ms]');
    %ylabel('V [mV] and I [nA]');
    %legend('Potential', 'Applied current');
    %replot

    %% export
    %print ("iaf_leaky-Be-potential.pdf", "-color", "-dpdf");
end;

%----------------------------------
% solve with Crank-Nicolson method
%----------------------------------
function ret = solveCN(lif)
    t = [];
    tlocal = 0.;
    v = lif.vleak;
    for i = 1 : lif.n-1
        % CN
        t(i)  = i * lif.dt;
        v = 1.0 / (1.0 + 0.5 * lif.dt / lif.taum) * (v + lif.dt/lif.taum *
            (0.5 * (lif.R * lif.I + lif.R * lif.I) + lif.vleak - 0.5 * v));

        % spike detected !
        if v >= lif.vpeak
            v  = lif.vreset; % membrane potential reset
            tlocal = t(i);
        end;
        
        % evaluate analytical solution
        vana = anaSol(lif, -tlocal + t(i));
        errV = v - vana;

        % save data in vectors for plotting
        ii(i)  = lif.I;
        vv(i)  = v;
        err(i) = errV;

        % print to standard cout
        %printf("%f\n", v);
    end;

    errV

    %% plot  solution
    %figure;
    %plot(t, vv, t, ii); grid;
    %xlabel('t [ms]');
    %ylabel('V [mV] and I [nA]');
    %legend('Potential', 'Applied current');
    %replot

    %% export
    %print ("iaf_leaky-CN-potential.pdf", "-color", "-dpdf");
end;

%--------------
% main program
%--------------
dt = 0.1;
for i = 1 : 5 
    lif.dt = lif.dt / 2;
    lif.n  = round(lif.T / lif.dt) + 1;
    %solveFE(lif);
    %solveBE(lif);
    solveCN(lif);
end;
