%%  Part 1:  Download and examine the IF model
%   To start, download the function:
%
%   IF.m
%
%   from http://math.bu.edu/people/mak/samsi/IF.m
%
%   Examine this function ... is it consistent with the model we proposed? 

%%  Part 2:  Execute the IF model
%   To run the IF model, we'll first need to define some parameters,

clear                   %Clear the workspace

I=0.1;                  %Set the parameter I = the injected current.
C=1;                    %Set the parameter C = the capacitance.
V0=0;                   %Set the initial voltage for the simulation.
Vth=2;                  %Set the voltage threshold.
Vreset=-1;              %Set the voltage reset.
N=10000;                %Set the # of steps to simulate.

%  Then,

[ V, t, n_spikes ] = IF( I, C, V0, Vth, Vreset, N );

plot(t,V)
xlabel('Time [s]');  ylabel('Voltage [mV]');

%IN LAB Q:  Adjust the parameter I.  What happens to V(t) if I=0?

%%  Part 3:  Determine the firing rate of the IF model.
%   Inspection of the simulation results suggests that the firing rate
%   changes with the injected current (parameter I).  Let's make a plot of
%   the firing rate versus the input current.  To do so, we'll simulate the
%   model for different values of I, keeping the other parameters fixed as
%   in Part 2.

I = (0:0.1:5);                 %Define a range of I values.
fr = zeros(size(I));            %... and a vector to hold spike counts.
for k=1:length(I)               %For each I, compute the # of spikes,
    [ V, t, n_spikes ] = IF( I(k), C, V0, Vth, Vreset, N );
    fr(k) = n_spikes / (t(end)); %... and save the firing rate.
end

plot(I,fr, 'r')
xlabel('I');  ylabel('Firing rate')

%IN LAB Q:  The resulting f-I curve appears to be a line.  What is the
%slope?

%% CHALLENGE QUESTION.
%
%  Implement and simulate the *leaky* integrate and fire model. Compute the
%  f-I curve for this model using the following parameter settings:
%
%  I=10;                   %Set the parameter I = the injected current.
%  C=1;                    %Set the parameter C = the capacitance
%  R=10;                   %Set the parameter R = the resistance
%  V0=-75;                 %Set the initial voltage for the simulation
%  Vth=20;                 %Set the voltage threshold.
%  Vreset=-90;             %Set the voltage reset.
%  Vrest=-70;              %Set the resting voltage
%  N=1000;                 %Set the # of stesp to simulate.

