% ********************************************************************
%   Quadratic Integrate-and-Fire model, solved numerically using Euler                                     
%   Math 430 * Math 635 * Lecture 7 * Sep 30, 2015 * Victor Matveev  
% ********************************************************************

clear;
b      =  1;    % Main parameter of quadratic I&F model: V' = V^2 - b + I
Vspike =  5;    % Spike peak = 5
Vreset = -5;    % Reset potential = -5
 
Iamp = 2;       % Current pulse amplitude
timeTotal = 70; % Total simulation time

% 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 > 10) .* (t < 60); 
f  = @(V, I) V.^2 - b + I;  % Rate of change of potential
                                           
h  = 0.1;               % Euler method time step
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) = -sqrt(b);        % Start the cell voltage at rest
n  = 1;                 % Initial value of array index

while n < length(T)
    n    = n + 1;
    V(n) = V(n-1) + h * f(V(n-1), I(n-1)); % Euler method 
    
    if V(n) > Vspike       % If threshold is reached, then...
        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)');
    title(['Current pulse amplitude: I=', num2str(Iamp)]);
    
    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)');
    title('Potential, V(t)');
    
    drawnow;
end