|
@@ -0,0 +1,144 @@
|
|
|
+function [STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
|
|
|
+%[STA_output]=BC_STA(loadData,filename,pixelSize,seed)
|
|
|
+%this function loads the voltage trace of the selected bipolar cell (BC)
|
|
|
+%and perorms a response-weighted average analysis analogous to the commonn
|
|
|
+%spike-triggered average (STA). For simplicity we refer to the analysis
|
|
|
+%STA,even though we work with a continous voltage signal and not spikes.
|
|
|
+%Inputs:
|
|
|
+% loadData = the folder path of the white noise files
|
|
|
+% filename = the filename of the BC of interest
|
|
|
+% pixelSize = the size of a square (checker) in monitor
|
|
|
+% pixels (see Documentation_Data.pdf)
|
|
|
+% seed = seed for the random number generator
|
|
|
+% Hz = stimulus update rate
|
|
|
+%
|
|
|
+%Outpus: STA_output = structure containing the STA, stimulus signal,
|
|
|
+% averaged membrane potential for calculating the output
|
|
|
+% nonlinearity
|
|
|
+%
|
|
|
+%
|
|
|
+% Note: the function is built only for cells recorded without frozen frames.
|
|
|
+% But additional comments are made for cells with frozen frames.
|
|
|
+% code by HS, last modified Feb 2021
|
|
|
+
|
|
|
+%% constants
|
|
|
+%OLED dimensions for regenerating the stimulus
|
|
|
+monitorXPixels=800; %x pixels of OLED monitor;
|
|
|
+monitorYPixels=600; %y pixels of OLED monitor;
|
|
|
+microPerPix=2.5; %one pixel=2.5 micrometer in size
|
|
|
+STA_TimeBack=1000; %time in ms you want to go back for generating the STA (to get smaller data I use here 1000ms, in the manuscript we use 2000ms).
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 1.load cell
|
|
|
+data=load([loadData,filename]); %load the file
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 2. regenerate the stimulus
|
|
|
+Nframes=length(data.ttls_ms); %number of frames shown (how many times the stimulus was updated)
|
|
|
+Nx=ceil(monitorXPixels/pixelSize); %number of squares shown X dimension
|
|
|
+Ny=ceil(monitorYPixels/pixelSize); %number of squares shown Y dimension
|
|
|
+nbrCheckers_total=Nx*Ny; %total number of checkers/square presented on the screen
|
|
|
+stimulus = recreateBinaryWhiteNoiseStimulus(Nx, Ny, Nframes, seed); %the function is available under: https://github.com/gollischlab/RecreateWhiteNoiseStimuliWithMatlab
|
|
|
+stimulus =stimulus(:,:,1:end-1);%take the last frame out
|
|
|
+%the output stimulus has 3 dimensions: Y,X (nbr of checkers on Y and X axis) and time. Easiest to understand the stimulus is by
|
|
|
+%plotting for example the first stimulus frame presented on the retina:
|
|
|
+%imagesc(stimulus(:,:,1));
|
|
|
+%Note: the stimulus values are drawn onto the
|
|
|
+%retina from the lower left corner, yet matlab displays it from the upper left corner.
|
|
|
+%Therefore the y axis has to be flipped when displaying.
|
|
|
+%set(gca,'YDir','normal');
|
|
|
+%For me the easiest is to directly flip the stimulus y dimension (and then I dont flip the axis when displaying):
|
|
|
+stimulus_origD = flipdim(stimulus,1);
|
|
|
+clear stimulus
|
|
|
+%Yet that might be personal preferences.Note: For generating the STA it does not
|
|
|
+%matter whether it is flipped or not. However if the STA is convolved with the natural movies it matters.
|
|
|
+%Further, if you have frozen and running frames, regenerate the stimulus separatly
|
|
|
+%for running and frozen frames (see Documentation_Data.pdf for the corresponding seeds).
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 3. make average per stimulus frame
|
|
|
+membranPotAVG=avgPerFrame(data.membPot,data.ttls_ms); % weights needed for computing the STA
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 4. Computing the STA
|
|
|
+%decide for numbers of frames you want to go back in time
|
|
|
+STA_fnbr=(ceil(STA_TimeBack/(1000/Hz))); % define the number of frame you want to go back in time
|
|
|
+%here I prefer to work with a 2D signal of the stimulus
|
|
|
+stimulus_2D = reshape(stimulus_origD, nbrCheckers_total, []);
|
|
|
+tic
|
|
|
+[~,STA_2D]=buildSTA ...
|
|
|
+ (stimulus_2D,membranPotAVG,STA_fnbr);
|
|
|
+toc
|
|
|
+%if you like the STA as a 3D e.g.:
|
|
|
+STA_3D=reshape(STA_2D,Ny,Nx, ...
|
|
|
+ size(STA_2D,2));% Y (monitor axis), X (monitor axis) and Z (STA_fnbr)
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 5. Make a zoom-in of the STA
|
|
|
+%cut a window of around 100 micrometer on each side from the maximum. The
|
|
|
+%zoom-in of the STA is better to avoid overfitting for the NL and
|
|
|
+%prediction.
|
|
|
+oneCheckInMicron=pixelSize*microPerPix;
|
|
|
+zoom_RF=ceil(100/oneCheckInMicron);
|
|
|
+
|
|
|
+%get the maximum of the STA to make a zoom-in
|
|
|
+[max_Color,indx_Zoom]=max(abs(STA_3D(:)));
|
|
|
+[posY,posX,posZ]=ind2sub(size(STA_3D),indx_Zoom);
|
|
|
+%cut a window around the STA ->make a zoom-in STA
|
|
|
+%zoom STA -> here you take the x-y coordinates
|
|
|
+rangeY_zoom = posY+(-zoom_RF:zoom_RF); rangeY_zoom = rangeY_zoom(rangeY_zoom>0 & rangeY_zoom<=size(STA_3D,1));
|
|
|
+rangeX_zoom = posX+(-zoom_RF:zoom_RF); rangeX_zoom = rangeX_zoom(rangeX_zoom>0 & rangeX_zoom<=size(STA_3D,2));
|
|
|
+STA_3D_zoom=STA_3D(rangeY_zoom,rangeX_zoom,:); %take the data from the 2000ms STA
|
|
|
+%reshape the STA
|
|
|
+STA2D_zoom=reshape(STA_3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(STA_3D_zoom,3));
|
|
|
+%get the zoom in of the stimulus signal
|
|
|
+stimulus3D_zoom=stimulus_origD(rangeY_zoom,rangeX_zoom,:); %take the data from the 2000ms STA
|
|
|
+stimulus2D_zoom=reshape(stimulus3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(stimulus3D_zoom,3));
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 5. visualize STA: plot the frame with maximal intensity value
|
|
|
+%I also like to plot all/multiple STA frames around the maximum.
|
|
|
+colorLimits = [-max_Color,max_Color]; %mainly relevant if you plot multiple frames and want to keep the same colorbar
|
|
|
+%plot
|
|
|
+figure;
|
|
|
+imagesc(STA_3D(:,:,posZ))
|
|
|
+caxis(colorLimits);
|
|
|
+colormap('gray'); %how the sitmulus was shown
|
|
|
+axis equal; axis tight;
|
|
|
+title(['STA best frame, filename: ' filename])
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 6. add variables into the output
|
|
|
+STA_output.STA_2D=STA_2D;
|
|
|
+STA_output.STA2D_zoom=STA2D_zoom;
|
|
|
+% STA_output.STA=STA; %not normalized
|
|
|
+STA_output.stimulus2D_zoom=stimulus2D_zoom;
|
|
|
+STA_output.Nx=Nx; %for reshapeing to 3D structure if needed
|
|
|
+STA_output.Ny=Ny;
|
|
|
+STA_output.membranPotAVG=membranPotAVG;
|
|
|
+STA_output.STA_fnbr=STA_fnbr;
|
|
|
+STA_output.filename=filename;
|
|
|
+STA_output.rangeX_zoom=rangeX_zoom;
|
|
|
+STA_output.rangeY_zoom=rangeY_zoom;
|
|
|
+STA_output.pixelSize=pixelSize;
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 7. additional comments:
|
|
|
+%From the zoom STA, I seperate the response into the highest-ranked spatial and temporal
|
|
|
+%components by singular-value decomposition (see Mansucript Schreyer & Gollisch, 2021)
|
|
|
+%example code for the SVD:
|
|
|
+%[sp,sigma,temp] = svd(STA,'econ'); %use the non-normalized STA here!
|
|
|
+%sp=spatial component and temp= the temporal component.
|
|
|
+%the highest-ranked spatial and temporal component has to be defined e.g. by manually looking
|
|
|
+%at the components.
|
|
|
+%From here, I fitted a 2D Gaussian to the best spatial component.
|
|
|
+%Important for nonlinearity calculation:
|
|
|
+%To avoid noise contributions from pixels outside the receptive field, I reduced
|
|
|
+%the number of elements in the STA by setting pixel values of the STA outside
|
|
|
+%the 3-sigma contour of the Gaussian fit to zero and re-separating the spatiotemporal
|
|
|
+%receptive field within this window into the highest-ranked spatial and temporal components
|
|
|
+%by singular-value decomposition. Each component was then normalized to unit Euclidean norm. E.g.
|
|
|
+%sp_bestnorm=normSTA(sp_best);
|
|
|
+%temp_bestnorm=normSTA(temp_best);
|
|
|
+
|
|
|
+end
|