% ********************************************************************
%   Integrate-and-Fire model, solved numerically using Euler method                                     
%   Math 430 * Math 635 * Lecture 8 * Oct 2, 2015 * Victor Matveev  
% ********************************************************************

clear;
VR     = -65;    % Resting potential = -65 mV
Vth    = -50;    % Threshold potential = -50 mV
Vspike =  10;    % Spike peak = 10mV (spike inserted "by hand" when V>=Vth)
Vreset = -75;    % Reset potential = -75 mV
 
tau  = 30;       % Membrane time constant in ms
R    = 4;        % Input resistance in MegaOhm (MOhm*nA = mV: convenient!)
Iamp = 5;        % Current pulse amplitude in nA
timeTotal = 700; % Total simulation time in ms

% Logical operations can be used to create a dis-continuous pulse function:  
% Current pulse of amplitude Iamp and width of 500ms as a function of time:

Current = @(t) Iamp * (t > 100) .* (t < 600);                                         
f = @(V, I) -(V - VR - R*I)/tau;  % Rate of change of potential

h  = 1;                 % Euler method time step, in ms
T  = 0 : h : timeTotal  % this is the time value array; step size = dt
I  = Current(T);        % Pre-compute current vs t, and store in array I
V(1) = VR;              % Start the cell voltage at rest
n  = 1;                 % Initial value of array index

while n < length(T)
    n    = n + 1;
    VRI  = VR + R * I(n);  % Equilibrium potential depends on the current
    V(n) = V(n-1) + h * f(V(n-1), I(n-1));  % Euler method 
    
    if V(n) > Vth       % If threshold is reached, then...
        V(n) = Vspike;  % ...insert a spike "by hand"
        n    = n + 1;   % ...advance by one time step, and...
        V(n) = Vreset;  % ...reset the potential to Vreset
    end
   
    subplot(2,1,1);    % Top figure panel: I(t), nA
    plot( T(1:n), I(1:n), 'r-', 'linewidth', 2); 
    axis([T(1) T(end) 0 Iamp]);  % Specify coordinate ranges
    xlabel('t, ms'); ylabel('I(t), nA');
    title(['500 ms current pulse, I=', num2str(Iamp), 'nA']);
    
    subplot(2,1,2);    % Bottom figure panel: V(t), mV
    plot( T(1:n), V(1:n), 'b-', 'linewidth', 2);
    axis([T(1) T(end) Vreset Vspike]); % Specify coordinate ranges
    xlabel('t, ms'); ylabel('V(t), mV');
    title('Potential, V(t) (mV)');
    
    drawnow;
end