123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 |
- function [STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
- %[STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
- %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 March 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 ->needed for doing the average per frames see avgPerFrame
- %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
- %For simplification and to keep the code simple, I cut a window of around
- %100 micrometer on each side from the maximum. A zoom-in of the STA is better
- %to avoid overfitting for the NL and prediction.See Methods part of the
- %manuscript and comment-section 7 below for how it is done in the manuscript.
- 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,:);
- %reshap to 2D
- 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 and prediction 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. To compute the nonlinearity, each component
- %was normalized to unit Euclidean norm. E.g.
- %sp_bestnorm=normSTA(sp_best);
- %temp_bestnorm=normSTA(temp_best);
- end
|