% I_and_F_slow_sodium
%
% Spring 2008

% Integrate-and-fire model with slow post-synaptic-spike-triggered Na+ conductance
%   (a simple form of "spike-response" model)

clear
epoch = 10;                         % length of epoch (sec)
bin_length = .0001;                 % bin length (sec)

E_Na = 0;                           % Na+ equilibrium potential (V)
delta_r_m_g_Na = .0091;             % increment of r_m*g_Na following spike emission (dimensionless)
tau_Na = .93;                       % decay time constant for r_m_g_Na (sec)

time_axis = 0:bin_length:epoch;     % time axis (bin units)
N = length(time_axis);              % number of bins in spike train

Vrest = -.065;                      % resting potential (V)
tau_m = .030;                       % membrane time constant (sec)
Vth = -.050;                        % firing threshold (V)
Vreset = -.065;                     % reset potential (V)
Rm = 9*10^7;                        % membrane resistance (Ohm)
Ie  = 10^-9;                        % electrode current (A)
window = .1/bin_length:1/bin_length;% current-injection window (bin units)

current = zeros(1,N);               % electrode current will be zero except in window
current(window) = Ie;               % set electrode current at specified value in window
V = Vrest*ones(1,N);
r_m_g_Na = zeros(1,N);

for n = 2:N
    if V(n-1) < Vth

        V_derivative = (Vrest - V(n-1) + r_m_g_Na(n-1)*(E_Na - V(n-1)) + Rm*current(n))/tau_m;
        V(n) = V(n-1) + bin_length*V_derivative;

        r_m_g_Na_derivative = - r_m_g_Na(n-1)/tau_Na;   
        r_m_g_Na(n) = r_m_g_Na(n-1) + bin_length*r_m_g_Na_derivative;

    else
        V(n-1) = .1;        % spike
        V(n) = Vreset;      % reset potential
        r_m_g_Na(n) = r_m_g_Na(n-1) + delta_r_m_g_Na;
    end
end


figure(1)
subplot(2,1,1)
plot(time_axis,V);
subplot(2,1,2)
plot(time_axis,r_m_g_Na);