%% iEIF - Inactivating Expoential Integrate and Fire Model 

% A reduced model of sodium channel inactivation and its impact on the 
% spike threshold

% References
%-------------------------------------------------------------------------%
% 1.  Platkiewicz J. and Brette, R. (2011) Impact of Fast Sodium Channel 
% Inactivation on Spike Threshold Dynamics and Synaptic Integration. PLoS
% Comput. Biol. 7(5): e1001129

% 2.  Platkiewicz, J. and Brette, R (2010) A Threshold Equation for Action 
% Potential Initiation. PLoS Comput. Biol. 6(7): e1000850 
%-------------------------------------------------------------------------%

% Written by S.E.C
% January 2018

%% Set Figure Properties

clear 
clf
close all
clc 

width = 2;      % Width in inches
height = 1.75;  % Height in inches
alw = 0.75;     % AxesLineWidth
fsz = 11;       % Fontsize
lw = 1;         % LineWidth
msz = 8;        % MarkerSize

set(0,'defaultLineLineWidth',lw)   % set the default line width to lw
set(0,'defaultLineMarkerSize',msz) % set the default line marker size to msz
set(0,'defaultLineLineWidth',lw)   % set the default line width to lw

% Set the default Size for display
defpos = get(0,'defaultFigurePosition');
set(0,'defaultFigurePosition', [defpos(1) defpos(2) width*100, height*100]);

% Set the defaults for saving/printing to a file
set(0,'defaultFigureInvertHardcopy','on')
set(0,'defaultFigurePaperUnits','inches')
defsize = get(gcf, 'PaperSize');
left = (defsize(1)- width)/2;
bottom = (defsize(2)- height)/2;
defsize = [left, bottom, width, height];
set(0, 'defaultFigurePaperPosition', defsize)
set(0,'DefaultAxesFontSize', 10)
set(0,'DefaultAxesLineWidth', 1)

%% Initialize variables and set model parameters

dt = 0.001;   % Time step (ms)     
tspan = 250;  % Input current duration (ms)

steps  = ceil(tspan/dt);
t = 0:dt:(steps-1)*dt;

gNa = 0.036;  % Sodium conductance (mS/cm^2)
ENa = 50;     % Sodium reversal potential (mV)
gL = 0.18;    % Leak conductance (mS/cm^2)
eL = -55;     % Leak current reversal potential (mV)
R = 1/gL;     % Membrane specific resistance (MOhm)
tau = 5;      % Membrane time constant(ms)
tau_hf = 15;  % Fast inactivation time constant (ms)
tau_hs = 500; % Slow inactivation time constant (ms)
V_r = -70;    % Reset voltage (mV) 
sig = 0.5;    % Gaussian noise amplitude (mV)
Va = -38.6;   % Sodium activation reversal potential (mV)
ka = 4;       % Na+ activation slope factor (mV)

% Compute V_tresh (mV)
V_thresh = Va - ka*log((gNa/gL)*(ENa-Va)/ka);
disp(V_thresh)

% Input current (nA)
p = 50;
amp = 3.8;
I = [zeros(p/dt,1); amp*ones((tspan-2*p)/dt,1); zeros(p/dt,1)];

%% Simulate iEIF model

% Initialize variables
hf = zeros(steps,1);
hf(1) = 1;
hs = zeros(steps,1);
hs(1) = 1;
V = zeros(steps,1);
V(1) = -70;
theta = zeros(steps,1);
binary = zeros(steps,1);

% Stochastic Euler Method
for i=1:steps-1       
     
    V(i+1) = V(i) + (dt/tau)*(hf(i)*hs(i)*exp((V(i)-V_thresh)/ka) + ...
                  + (eL-V(i)) + R*I(i)) + sig*sqrt(dt)*randn;                 
    hf(i+1) = hf(i) + (dt/tau_hf)*(h_inf(V(i)) - hf(i));
    hs(i+1) = hs(i) + (dt/tau_hs)*(h_inf(V(i)) - hs(i));
    theta(i+1) = V_thresh - ka*(log(hs(i)) + log(hf(i))); 
    
    if V(i+1) >= theta(i+1)
        binary(i+1) = 1; 
        V(i+1) = V_r;
    end    
      
end

thresholds = theta(binary>0);
isi = diff(find(binary)*dt);
mu_isi = mean(isi);
disp([num2str(1/mu_isi*1000),' Hz'])

%% Display the model results

fig = figure(1);
%-------------------------------------------------------------------------%
subplot(5,1,1:2)
plot(t,V,'k')
hold on
set(gca,'Ylim',[V_r 15])
hold on 
vline(find(binary)*dt,'k')
set(gca,'Xlim',[0 tspan],'Xtick',[],'Ylim',[-80 20],...
   'Ytick',[-80 20],'TickLength',[0 0])
box on
%-------------------------------------------------------------------------%
theta2 = ones(1,length(t))*theta(find(binary,1,'last'));
theta1 = theta(find(binary,1));
subplot(5,1,4:5) 
h = area(theta2,theta1,'FaceColor',[0.9 0.9 0.9]);
h.BaseLine.Color = [0.9 0.9 0.9];
txt = annotation('textbox', [0.9, 0.31, 0, 0], 'string','4');
txt.Color = [0.7 0.7 0.7]; 
hold on
plot(t,theta,'k')
set(gca,'Xlim',[0 tspan],'Xtick',[0 tspan],'Ylim',[-45 -20],...
'Ytick',[-45 -20],'TickLength',[0 0])  
box on
%-------------------------------------------------------------------------%
subplot(5,1,3)
plot(t,I,'k')
set(gca,'Xlim',[0 tspan],'Xtick',[],'Ylim',[0 round(amp)+1],'TickLength',[0 0])
set(gca,'Visible','off')
B = get(gca,'position');        % B is a 4D vector [lef bot wid hei]
B(1,4) = B(1,4)/2;              % reduce the height by half
B(1,2) = B(1,2) + 1.1*B(1,4);   % change the vertical position
set(gca,'position',B)     


figure(2)
subplot(2,1,1)
plot(V,hf,'k')
set(gca,'Xlim',[-75 -30],'Xtick',[-80 -30],'Ylim',[-0.05 1.05],'Ytick',[0 1],'TickLength',[0 0])
subplot(2,1,2)
plot(V,hs,'k')
set(gca,'Xlim',[-75 -30],'Xtick',[-80 -30],'Ylim',[0.6 1.05],'Ytick',[0.6 1],'TickLength',[0 0])

figure(3)
for k = 1:length(thresholds)
    clear temp
    temp = zeros(length(thresholds)-k,1);
    for j = 1:length(thresholds)-k
    temp(j) = thresholds(end-(j-1))-thresholds(end-(j-1)-k);
    end
    plot(k,mean(temp),'ko')
    hold on
end
set(gca,'Xlim',[1 15],'TickLength',[0 0])

%% Boltzman Na inactivation function

function aH = h_inf(V)
Vi = -63;     % Sodium inactivation reversal potential (mV)
ki = 6;       % Inactivation slope
aH = 1./(1+exp((V-Vi)/ki));
end

