%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation of (leaky) integrate-and-fire neuron
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all

%% Model Parameters
tau=15;                % time constant (ms)
AP_threshold=-30;      % AP threshold (mV)
I_ext=40.05;           % constant external input
noise_sd=0;            % number of SDs for gaussian input noise (mV)
rest_potential=-70;    % resting state potential (mV)
tspan=[1,750];        % integration interval
num_steps_per_ms=1;    % num of steps per ms

%% Set important variables
curr_voltage=rest_potential;  % initial membrane voltage (mV)
tau=tau*num_steps_per_ms; %update tau constant for step size
T=tspan(1):(1/num_steps_per_ms):tspan(2); %time vector (ms)

AP=zeros(length(T)); %initialize spike vector
voltage=zeros(length(T)); %initialize voltage vector
counter=0; %loop counter

%% Integration with Euler method
for t=T;
    counter=counter+1; %advance counter
    I_ext_curr=I_ext +randn*noise_sd; %noisy (gaussian) input
    
    if curr_voltage>AP_threshold %Does the neuron currently fire?
        curr_voltage=rest_potential; %reset voltage
        AP(counter)=1;  %mark AP at current time
    
    else %not firing - update voltage:
        curr_voltage=curr_voltage+... %the current voltage
            (rest_potential-curr_voltage)./tau+... %minus leak
            I_ext_curr/tau; %plus (noisy?) input
    end
    
    voltage(counter)=curr_voltage; %store current voltage
end

%% Plot
%lower subplot with voltages
subplot('position',[0.13 0.13 1-0.26 0.6])
plot(T,voltage);
hold on; 
plot(tspan,[AP_threshold AP_threshold],'--');
xlabel('time [\tau]');
ylabel('u(t)')

%upper subplot with spikes
subplot('position',[0.13 0.8 1-0.26 0.1])
plot(T,AP,'.','markersize',20);
axis([tspan(1), tspan(2), 0.5, 1.5])
set(gca,'xtick',[],'ytick',[])
ylabel('Action Potentials')
hold off