|
@@ -1,144 +0,0 @@
|
|
|
-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
|