%% IF neurons, multiple

% This script simulates a number of integrate and fire neurons, each of 
% which receives a constant input, plus a noise component common to all 
% neurons (representing the global network state, for example), plus an 
% independent noise component for each neuron (reflecting its unique input  
% connections coming from other neurons)
% 
% The neurons are characterized by a resting membrane potential, a spike
% generation threshold, a reset value for the membrane potential after a
% a spike, and the membrane time constant. The input resistance is set to 1
% for convenience.
%
% In the end, the input to the cells is plotted as well as the resulting
% raster plot with the spikes that were elicited for each neuron.
%
% Based on code by Mark van Rossum
% Matthias Hennig and Theoklitos Amvrosiadis, 2018



%% Model setup

close all % Closes all open figures

% Basic variables of the model
ncells = 50; % Number of cells
vrest = -70; % Resting membrane potential (mV)
vreset = -70; % Reset voltage
vthr = -50; % Spike threshold
taum = 20; % Membrane time constant (ms)

sdcommon = 100; % Noise, common to all cells
sd = 50; % Independent noise for each neuron

dt = 0.1; % Time-step duration (ms)
endtime = 1e3; % End of simulation (ms)
nsteps = endtime/dt; % Number of timesteps

% Create vectors to store the data
spikes = zeros(ncells,nsteps); % Spikes, binary
vmem = vrest*ones(ncells,nsteps); % Membrane voltage (for each timestep)
stimhist = zeros(1,nsteps); % Stores input current for plotting

vmem(:,1) = vrest; % Initialization value


%% Simulation run

for itime = 2:nsteps % Iterate for every timestep until the end of the simulation
    
    istim = 20; % Stimulus(input) amplitude (nA)
    
    % Adjust the input amplitude halfway through the simulation
    if (itime > nsteps/2) 
        istim = 40;
    end
    
	istim = istim + sdcommon*randn; % Add common noise for all cells
	stimhist(itime)= istim;
    itot = istim*ones(1,ncells) + sd*randn(1,ncells); % Total input for each cell, including independent noise
    vmem(:,itime) = vmem(:,itime-1) + dt/taum*(-(vmem(:,itime-1)-vrest) + itot'); % Calculate membrane voltage for current timestep

    spx = find(vmem(:,itime)> vthr); % Find which cells, if any, have reached spike threshold      
    spikes(spx,itime) = 1; % Store spikes
    vmem(spx,itime) = vreset; % Reset membrane voltage to reset value
    
end

fr = mean(sum(spikes,2)/(endtime/1e3)); % Average firing rate(Hz) across all neurons 

%% Plotting the results

% First, the raster plot of the spikes
for icell = 1:ncells
    
    spx = find(spikes(icell,:)); % Get the spike times for each cell
    spy = icell;

    for isp = 1:length(spx)
        subplot(4,1,[1,3])
        plot([spx(isp), spx(isp)], [icell-0.5, icell+0.5], 'k')
        hold on
    end
end
hold off

set(gca,'xtick',[])
ylabel('Neuron no.', 'FontSize', 10)
ylim([0, ncells+1])
title('Spikes')

% Second, the input current to the cells
stimhist = smooth(stimhist, 80)/10; % (rescale as well)
subplot(4,1,4)
plot((1:nsteps)*dt, stimhist, 'r')
ylim([0.8*min(stimhist(2:end-1)), 1.2*max(stimhist(2:end-1))])
yticks(floor(0.8*min(stimhist(2:end-1))):3:floor(1.2*max(stimhist(2:end-1))))
xlabel('Time(ms)', 'FontSize', 10)
xt = get(gca, 'XTick');
set(gca, 'FontSize', 10)
ylabel('Input current(10x nA)', 'FontSize', 10)

title('Input amplitude')

% Finally, plot the ISI histogram
figure(2)
isi = []; 

for icell = 1:ncells
    spx = find(spikes(icell,:));
    cellisi = diff(spx);
    isi = [isi cellisi];
end

isi = isi*dt;

histogram(isi, 'Normalization', 'probability')
title('ISI distribution')
xlabel('ISI (ms)')
ylabel('fraction/bin')
