%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrate and Fire Model
% Odelia Schwartz
% 2016

% Simple version of Integrate and Fire model. 
% dV/dt = - V/RC + I/C
% Voltage is determined here by directly solving the differential equation.

clear all;

deltaT=1; 				% 1 msec time steps
tend = 100;

% Define the range of times for the simulation.
times=[1:tend]; 			% 100 msec simulation

% Initialize membrane potential at rest, no spikes
Vstart = 0;
spikes=zeros(size(times));

% Loop through time, updating the membrane potential and adding
% spikes whenever the membrane potential crosses the spike
% threshold.

I = 10;   % current
R = 50;   % resistance
C = 6;    % capacitance
k = R*C;  % time constant
V(1) = Vstart;
Vthresh = 25; % threshold for spiking
t = 0;

for i=1:tend-1

  % update membrane potential
  V(i+1) = Vstart*exp(-times(t+1)./k) + R .* I .* (1 - exp(-times(t+1)./k));
  t = t+1;

  % If above spike threshold, generate a spike and reset the
  % membrane potential to Vleak
  if V(i+1) > Vthresh 
    spikes(i+1) = 1;
    V(i+1) = Vstart;
    t=0;
  end

end

% Plot membrane potential and spike train
subplot(2,1,1)
plot(times,V)
xlabel('Time (msec)')
ylabel('Membrane potential (mV)')
subplot(2,1,2)
stem(times,spikes, '-o')
xlabel('Time (msec)')
ylabel('Neural response')

%% To do
% Try changing the parameters and re-running:
% Increase I
% Set I to 0 and Vstart to 20
% Decrease R
% Change C

