InputCurrents.m

clear all
close all

%% Figure dimensions

figure(1) 
    set(gcf,'PaperUnits','centimeters') %Setting the units of the figure on paper to centimeters.
    xSize = 13; 
    ySize = 13; 
    %Size of the figure 

    xLeft = (21-xSize)/2;
    yTop = (30-ySize)/2; 
    %Coordinates to center the figure on A4-paper
    set(gcf,'PaperPosition',[xLeft yTop xSize ySize])
    %This command sets the position and size of the figure on the paper to the desired values. 
    set(gcf,'Position',[0.5 0.5 xSize*50 ySize*50]) 
    set(gcf, 'Color', 'w');

%% Parameter Initialization:

duration=5000;  % duration in ms
dt=0.1;         % simulation time step.
tau=50;         % Filter time for the input. 
tRef=5;         % Refractory period for the spike trains. 
nn=100;         % Number of spiketrains we seek to create later. 
spikebin=5;     % Number of ms per PSTH bin.

Timevector=(0.1:dt:duration);               % A vector of time in ms, in steps of dt.
WhiteNoise=rand(1,length(Timevector))-0.5;  % uniform white noise drawn from +/- 0.5 as long as the time vector.
FilteredWhiteNoise=WhiteNoise.*0;           % an empty vector which we will use to create the time-filtered input.
SpikeTrains=zeros(nn,length(Timevector));   %A Matrix that will hold all spiketrains.
PloTrains=SpikeTrains;                      % This is just a plotting variable to oercome a screen resolution problem in matlab.

avetrain=0;     % A counter to calculate the average firing rate.
tslt=0;         % (== t(ime)s(ince)l(ast)(t)oggle (this serves as a Boolean for the sparsification of the input signal.

tsls=zeros(nn,1);                               % (== t(ime)s(ince)l(ast)(s)pike (to keep track of the refractory period of each spike train)
BinnedSpikeTrains=zeros(1,duration/spikebin);   % a vector to create a PSTH with binwidth ?spikebin? from the spike trains.

%% Making the time-filtered white noise signal:

for t=2:duration/dt FilteredWhiteNoise(t) = WhiteNoise (t) - ...
            (WhiteNoise (t) - FilteredWhiteNoise(t-1))*exp(-dt/tau);
end

%% This routine changes the signal trace ?FilteredWhiteNoise? by a ?exp(-dt/tau)? fraction of the difference between the signal and a random number.

FilteredWhiteNoise=FilteredWhiteNoise./max(FilteredWhiteNoise);
%Normalize to a maximum value of 1.

%% Plotting:

figure(1) 
subplot(4,1,1) 
plot(Timevector, FilteredWhiteNoise) 
axis([0 duration -1 1]) 
x=sprintf('Time Filtered White Noise (FWN)'); 
title (x)

%% Normalize and Rectify:

FilteredWhiteNoise=FilteredWhiteNoise.*(500*dt/1000); % Normalizes the trace to a peak value of 500Hz*dt (=0.05).
FilteredWhiteNoise(FilteredWhiteNoise<0)=0; %Sets all negative values of ?FilteredWhiteNoise? to 0.

%% Plotting:

subplot(4,1,2) 
    hold on
    plot(Timevector,FilteredWhiteNoise, 'b', 'LineWidth', 1.1)

%% Sparsifieing the Input Signals:

% This routine goes through the signal trace and deletes entire ?activity bumps? if certain conditions are fullfilled: 
toggle = 0; 
tslt=0;

for d=1:duration/dt-1
    
    % Routine becomes active (sets toggle == 1) if the signal is ==0, and the toggle is off (==0) and has been off for at least 1 ms:
    if(FilteredWhiteNoise(d)==0 && toggle==0 && (d-tslt>10))
        toggle=1; % toggle set 
        tslt=d; % ?refractory? for toggle is set
    end
    
    % If routine active, every signal value is set to zero:
    
    if (toggle==1)
        FilteredWhiteNoise(d) = 0; % If routine has been active for longer than 0.5 ms, and the signal is 0, routine becomes inactive:    
        if (FilteredWhiteNoise(d+1)==0 && (d-tslt>5)) 
            toggle=0;
        end
    end
end

%% Plotting:
subplot(4,1,2) 
    hold on 
    plot(Timevector, FilteredWhiteNoise, 'r')
    axis([0 duration -0.005 0.05]) 
    title ('Rectified & callibrated (blue) and sparsened (red) FWN')

%% Adding background firing rate:

FilteredWhiteNoise=FilteredWhiteNoise+(5*dt/1000);

% This is adjusted so that without any FilteredWhiteNoise the firing rate is 5 Hz*dt (0.0005).

%% Creating 100 spike trains:
for i=1:nn 
    for t=1:duration/dt
        if (tsls (i) <= 0) % Allows potential spike if refractory period has subsided
            if(rand<FilteredWhiteNoise(t))
                SpikeTrains (i,t) = 1;%Firesifrandomvariable<?FilteredWhiteNoise?. 
                tsls (i) = tRef;% Sets the absolute refractory period.
                avetrain=avetrain+1;% Counts the total number of spikes.
                
                if(duration/dt-t>25)% This is just a plotting routine. 
                    PloTrains (i,t:t+25)=1;%(Spike is elongated for plotting.)
                end
            end
        else
            tsls (i)=tsls (i) -dt;% subtracts dt from refractory counter if it is still >0
        end
    end
end
avetrain=avetrain/(nn*duration/1000); %Calculates the average firing rate in Hz.
         

%% Plotting:
subplot(4,1,3) 
    imagesc(-PloTrains) 
    colormap(gray) 
    title ('100 Spiketrains')

%% Recording a PSTH / Binned Input rate:

% This bins the spikes into bins and calculates the instantaneous firing rate in Hz.
for i=1:(duration/spikebin)-1 
    BinnedSpikeTrains(i)=...
        sum(sum(SpikeTrains(:,((i-1)*(spikebin/dt))+1:(i*(spikebin/dt)))));
end

BinnedSpikeTrains= (BinnedSpikeTrains*(1000/spikebin))/nn;

%% Plotting:
subplot(4,1,4)
    plot(BinnedSpikeTrains)
    x=sprintf('Average input rate for 1 excitatory channel, \%3.2f Hz, peak \%3.2f Hz', avetrain, max(BinnedSpikeTrains));
    title (x)



%%End of code.