|
|
@@ -0,0 +1,384 @@
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+function Analyze_Spatial_Contrast_Adaptation_Stimulus(varargin)
|
|
|
+%
|
|
|
+%
|
|
|
+%%% Analyze_Spatial_Contrast_Adaptation_Stimulus %%%
|
|
|
+%
|
|
|
+% This function analyzes the response of the ganglion cells to the spatial
|
|
|
+% contrast adaptation stimulus. This is a supplement code for the Khani & Gollisch
|
|
|
+% (2017) study about the diversity of spatial contrast adaptation in the mouse retina.
|
|
|
+% This function accompanys dataset of Khani and Gollisch (2017). The dataset is available at:
|
|
|
+% https://gin.g-node.org/gollischlab/Khani_and_Gollisch_2017_RGC_spiketrains_for_spatial_contrast_adaptation
|
|
|
+% This function first loads the data in a format that is provided in by the
|
|
|
+% supplementary database of the study. It then generate PSTHs, spike-triggered
|
|
|
+% average (STA) and nonlinearity functions for the selected cell.
|
|
|
+%
|
|
|
+% ===============================Inputs====================================
|
|
|
+%
|
|
|
+% datapath : path to the folder in the database.
|
|
|
+% stimulusnumber : number to select correct stimulus.
|
|
|
+% cellnumber : an example cell number from the list of good cells.
|
|
|
+% clusternumber : a cluster number for the selected cell.
|
|
|
+%
|
|
|
+%================================Output====================================
|
|
|
+%
|
|
|
+% PSTH, STAs and nonlinearities for the selected cell.
|
|
|
+%
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+%---------- main-function ----------%
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+% first load the data for an example cell
|
|
|
+[dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin{:});
|
|
|
+% analyze the data
|
|
|
+ca_data = sca_anylsis(dp, stimulusnumber, cellnumber, clusternumber);
|
|
|
+% plotting
|
|
|
+plot_sca(ca_data, cellnumber, clusternumber);
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+%---------- sub-functions ----------%
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function expdata = sca_anylsis(dp, stimulusnum, cellnum, clusternum,varargin)
|
|
|
+%
|
|
|
+expdata = load_data_for_selected_cell(dp, stimulusnum, cellnum, clusternum);
|
|
|
+% set all the required parameters
|
|
|
+para = expdata.stimpara;
|
|
|
+para.refreshrate = 60;
|
|
|
+para.upsampleratio = ceil(1e3/para.refreshrate);
|
|
|
+para.numNLbins = 40;
|
|
|
+para.meanIntensity = 0.5;
|
|
|
+para.statime = fliplr(linspace(0,(1e3/60*para.Nblinks)*(para.numNLbins/para.upsampleratio),para.numNLbins));
|
|
|
+Nruns = floor((length(expdata.frametimes)-1)/(para.Nblinks*para.switchFrames));
|
|
|
+para.Nframes = Nruns*para.Nblinks*para.switchFrames;
|
|
|
+para.Nrepeats = floor(Nruns/2); % number of runs per trial;
|
|
|
+para.binlength = 250 / 1e3; % 250 ms each bin
|
|
|
+% for the PSTH
|
|
|
+para.low_contduration = para.switchFrames ./ para.refreshrate;
|
|
|
+para.high_contduration = para.switchFrames ./ para.refreshrate;
|
|
|
+para.psthtime_low = linspace(0, para.low_contduration ,(para.low_contduration / para.binlength)+1);
|
|
|
+para.psthtime_high = linspace(0, para.high_contduration ,(para.high_contduration / para.binlength)+1);
|
|
|
+para.trialduration = para.low_contduration + para.high_contduration;
|
|
|
+para.psthtimelh = linspace(0, para.trialduration,(para.trialduration / para.binlength));
|
|
|
+
|
|
|
+% loading clusters and frametimings
|
|
|
+ft = expdata.frametimes; % work in seconds
|
|
|
+ftimes = ft(1:para.Nframes);
|
|
|
+% get frame times per trial
|
|
|
+ftimes = reshape(ftimes,para.switchFrames,Nruns)';
|
|
|
+ftimes_high = ftimes(2:2:end,:); ftimes_high = ftimes_high(1:para.Nrepeats,:);
|
|
|
+ftimes_low = ftimes(1:2:end,:); ftimes_low = ftimes_low(1:para.Nrepeats,:);
|
|
|
+
|
|
|
+% getting the stimulus right!
|
|
|
+stimulus = Spatial_Contrast_Adaptation_Stimulus(para.seed,para.Nframes,para.switchFrames,Nruns, para.meanIntensity);
|
|
|
+
|
|
|
+% binning spikes
|
|
|
+spkbins = (histc(expdata.spiketimes,ft(1:para.Nframes))); %#ok
|
|
|
+spkbins = reshape(spkbins,para.switchFrames,Nruns)';
|
|
|
+spkbins_high = spkbins(2:2:end,:); spkbins_high = spkbins_high(1:para.Nrepeats,:);
|
|
|
+spkbins_low = spkbins(1:2:end,:); spkbins_low = spkbins_low(1:para.Nrepeats,:);
|
|
|
+
|
|
|
+% setting output
|
|
|
+expdata.ftimes.high = ftimes_high;
|
|
|
+expdata.ftimes.low = ftimes_low;
|
|
|
+expdata.stimulus = stimulus;
|
|
|
+expdata.stimpara = para;
|
|
|
+
|
|
|
+% calculating rasters and psth
|
|
|
+expdata = calc_rasters_psth(expdata);
|
|
|
+
|
|
|
+% calculating STAs and nonlinearities for each contrast and location
|
|
|
+[locX_high.sta,locX_high.nlx,locX_high.nly,locX_high.time] = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locX_high, para);
|
|
|
+[locX_low.sta,locX_low.nlx,locX_low.nly,locX_low.time] = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locX_low, para);
|
|
|
+[locY_high.sta,locY_high.nlx,locY_high.nly,locY_high.time] = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locY_high, para);
|
|
|
+[locY_low.sta,locY_low.nlx,locY_low.nly,locY_low.time] = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locY_low, para);
|
|
|
+
|
|
|
+% setting output
|
|
|
+expdata.locXhigh = locX_high;
|
|
|
+expdata.locX_low = locX_low;
|
|
|
+expdata.locY_high = locY_high;
|
|
|
+expdata.locY_low = locY_low;
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function expdata = calc_rasters_psth(expdata)
|
|
|
+
|
|
|
+fthigh = expdata.ftimes.high;
|
|
|
+ftlow = expdata.ftimes.low;
|
|
|
+spks = expdata.spiketimes;
|
|
|
+para = expdata.stimpara;
|
|
|
+[ras_high, ras_low] = deal(cell(para.Nrepeats,1));
|
|
|
+delay = 50/1e3; % ignore first spikes after contrast switch
|
|
|
+
|
|
|
+for ii = 1:para.Nrepeats
|
|
|
+
|
|
|
+ ras_high{ii} = spks(and(spks >= fthigh(ii,1),spks <= fthigh(ii,end))) - fthigh(ii,1);
|
|
|
+ ras_low{ii} = spks(and(spks >= (ftlow(ii,1)+delay),spks <= ftlow(ii,end))) - (ftlow(ii,1)+delay);
|
|
|
+
|
|
|
+end
|
|
|
+rasters_high = CelltoMatUE(ras_high);
|
|
|
+rasters_low = CelltoMatUE(ras_low);
|
|
|
+
|
|
|
+psth_high = mean(histc(rasters_high',para.psthtime_high),2)*(1/para.binlength); %#ok
|
|
|
+psth_high = psth_high(1:end-1); psth_high(end) = psth_high(end-1);
|
|
|
+para.psthtime_high = para.psthtime_high(1:end-1);
|
|
|
+psth_low = mean(histc(rasters_low',para.psthtime_low),2)*(1/para.binlength); %#ok
|
|
|
+psth_low = psth_low(1:end-1); psth_low(end) = psth_low(end-1);
|
|
|
+para.psthtime_low = para.psthtime_low(1:end-1);
|
|
|
+
|
|
|
+expdata.rasters_high = rasters_high;
|
|
|
+expdata.rasters_low = rasters_low;
|
|
|
+expdata.psth_high = psth_high';
|
|
|
+expdata.psth_low = psth_low';
|
|
|
+expdata.psthlh = [psth_low; psth_high]';
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+function [staNorm,nlx,nly,timeVec] = calculate_sta_nonlinearity(runningbin, stimulus, para)
|
|
|
+
|
|
|
+toHzratio = para.Nblinks/para.refreshrate; % normalize the response to Hz
|
|
|
+
|
|
|
+rnstimEn = zeros(para.Nrepeats,para.numNLbins, para.switchFrames-para.numNLbins+1);
|
|
|
+for itrial = 1 : para.Nrepeats
|
|
|
+ stimtrial = stimulus(itrial,:);
|
|
|
+ hankelMat = hankel(1:para.numNLbins,para.numNLbins:para.switchFrames);
|
|
|
+ rnstimEn(itrial,:,:) = stimtrial(hankelMat);
|
|
|
+end
|
|
|
+runningbin = reshape(permute(runningbin(:,para.numNLbins:end,:),[1 3 2]),1,[]);
|
|
|
+rnstimEn = reshape(permute(rnstimEn, [2 1 3]), para.numNLbins,[]);
|
|
|
+
|
|
|
+sta = bsxfun(@rdivide,runningbin*rnstimEn',sum(runningbin,2));
|
|
|
+staNorm = bsxfun(@times,sta,1./sqrt(sum(sta.^2,2)));
|
|
|
+
|
|
|
+
|
|
|
+timeVec = (-(para.numNLbins-1):1:0)*toHzratio; %in seconds
|
|
|
+
|
|
|
+runningGenerators = staNorm * rnstimEn;
|
|
|
+
|
|
|
+[nly,nlx] = calculate_nonlinearity(runningGenerators', runningbin', para.numNLbins, toHzratio);
|
|
|
+
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+
|
|
|
+function [ values, centers,sevalues] = calculate_nonlinearity(linearOutput, spikes, nbins, dt)
|
|
|
+
|
|
|
+qtStep = 1/nbins;
|
|
|
+qtEdges = 0:qtStep:1;
|
|
|
+genEdges = quantile(linearOutput, qtEdges);
|
|
|
+
|
|
|
+if all(isnan(genEdges)) % this is for cases without spikes where everything is NaNs
|
|
|
+ values = zeros(size(genEdges));
|
|
|
+ sevalues = zeros(size(genEdges));
|
|
|
+else
|
|
|
+ [a, spikebins] = histc(linearOutput,genEdges); %#ok
|
|
|
+
|
|
|
+ values = accumarray(spikebins, spikes', [nbins+1 1], @sum);
|
|
|
+ sdvalues = accumarray(spikebins, spikes', [nbins+1 1], @std);
|
|
|
+
|
|
|
+ values = values./a/dt;
|
|
|
+ sevalues = sdvalues./sqrt(a)/dt;
|
|
|
+
|
|
|
+end
|
|
|
+values(end) = [];
|
|
|
+sevalues(end) = [];
|
|
|
+centers = genEdges(1:nbins)+diff(genEdges)/2;
|
|
|
+values = transpose(values);
|
|
|
+sevalues = transpose(sevalues);
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function [dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin)
|
|
|
+% first parse the user inputs
|
|
|
+p = inputParser(); % check the user options.
|
|
|
+p.addParameter('datapath', [], @(x) ischar(x) || isstring(x));
|
|
|
+p.addParameter('stimulusnumber', [], @isnumeric);
|
|
|
+p.addParameter('cellnumber', [], @isnumeric);
|
|
|
+p.addParameter('clusternumber', [], @isnumeric);
|
|
|
+p.addParameter('help', false, @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
|
|
|
+
|
|
|
+p.parse(varargin{:});
|
|
|
+% defualt parameters
|
|
|
+defpara = p.Results;
|
|
|
+
|
|
|
+if isempty(defpara.datapath)
|
|
|
+ dp = uigetdir('Select an experiment folder:');
|
|
|
+else
|
|
|
+ dp = defpara.datapath;
|
|
|
+end
|
|
|
+
|
|
|
+if isempty(defpara.cellnumber) || isempty(defpara.clusternumber)
|
|
|
+ listofgoodcells = load([dp,filesep,'list_of_good_cells.txt'],'-ascii');
|
|
|
+ cellabels = sprintfc('cell %g, cluster %g',listofgoodcells);
|
|
|
+ if numel(cellabels) > 1
|
|
|
+ idx = listdlg('PromptString',{'Select an example ganglion cell:'},...
|
|
|
+ 'SelectionMode','single','ListString',cellabels,'listsize',[250 250]);
|
|
|
+ else
|
|
|
+ idx = 1;
|
|
|
+ end
|
|
|
+ cellnumber = listofgoodcells(idx,1);
|
|
|
+ clusternumber = listofgoodcells(idx,2);
|
|
|
+else
|
|
|
+ cellnumber = defpara.cellnumber;
|
|
|
+ clusternumber = defpara.clusternumber;
|
|
|
+end
|
|
|
+
|
|
|
+if isempty(defpara.stimulusnumber)
|
|
|
+ stimnames = importdata([dp,filesep,'stimuli_names.txt']);
|
|
|
+ canames = stimnames(~(cellfun('isempty',(strfind(stimnames,'Contrast_Adaptation')))));
|
|
|
+
|
|
|
+ if numel(canames) > 1
|
|
|
+ idx = listdlg('PromptString',{'Select an example dataset:'},...
|
|
|
+ 'SelectionMode','single','ListString',canames,'listsize',[450 100]);
|
|
|
+ canames = canames{idx};
|
|
|
+ end
|
|
|
+ stimulusnumber = str2double(extractBefore(canames,'_'));
|
|
|
+else
|
|
|
+ stimulusnumber = defpara.stimulusnumber;
|
|
|
+end
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function expdata = load_data_for_selected_cell(dp, stimnum, cellnum, clusternum)
|
|
|
+%
|
|
|
+ftfolder = [dp,filesep,'frametimes'];
|
|
|
+if ~exist(ftfolder,'dir'), ftfolder = dp; end
|
|
|
+rasfolder = [dp,filesep,'spiketimes'];
|
|
|
+if ~exist(rasfolder,'dir'), rasfolder = dp; end
|
|
|
+paramfolder = [dp,filesep,'stimulusparameters'];
|
|
|
+if ~exist(paramfolder,'dir'), paramfolder = dp; end
|
|
|
+
|
|
|
+ftnames = dir([ftfolder, filesep, num2str(stimnum,'%02g'),'*_frametimings.txt']);
|
|
|
+rastername = [rasfolder, filesep, sprintf('%d_SP_C%d%02d.txt', stimnum, cellnum, clusternum)];
|
|
|
+
|
|
|
+warning('off','all'); % this is supress the warning about long names, hopefully nothing happen in between these lines!
|
|
|
+expdata.spiketimes = load(rastername, '-ascii');
|
|
|
+expdata.frametimes = load([ftfolder,filesep,ftnames.name], '-ascii');
|
|
|
+
|
|
|
+stimpath = dir([paramfolder,filesep,num2str(stimnum,'%02g'),'*_parameters.txt']);
|
|
|
+
|
|
|
+% now reading the txt file with parameters
|
|
|
+fid = fopen([paramfolder,filesep,stimpath.name]);
|
|
|
+tline = fgetl(fid);
|
|
|
+while ischar(tline)
|
|
|
+ tline = fgetl(fid);
|
|
|
+ if tline == -1, break; end
|
|
|
+ fn = extractBefore(tline,' = ');
|
|
|
+ if isempty(fn), continue; end
|
|
|
+ val = extractAfter(tline,' = ');
|
|
|
+ % to convert to double
|
|
|
+ if ~isnan(str2double(val))
|
|
|
+ val = str2double(val);
|
|
|
+ end
|
|
|
+ % to convert to logical
|
|
|
+ if strcmp(val,'true'), val = true; end
|
|
|
+ if strcmp(val,'false'), val = false; end
|
|
|
+ % store the values in a structure
|
|
|
+ stimpara.(fn) = val;
|
|
|
+end
|
|
|
+fclose(fid);
|
|
|
+
|
|
|
+expdata.stimpara = stimpara;
|
|
|
+expdata.name = extractBefore(ftnames.name,'_frametimings.txt');
|
|
|
+expdata.paraname = extractBefore(stimpath.name,'_parameters.txt');
|
|
|
+if strcmp(expdata.name , expdata.paraname)
|
|
|
+ expdata = rmfield(expdata,'paraname');
|
|
|
+end
|
|
|
+warning('on','all');
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function OutputMat = CelltoMatUE (input_cell)
|
|
|
+%
|
|
|
+inptsize = size(input_cell);
|
|
|
+maxLength = max(cellfun('length',input_cell));
|
|
|
+if inptsize(2) > inptsize(1)
|
|
|
+ OutputMat = cell2mat(cellfun(@(x)cat(1,x,nan(maxLength-length(x),1)),input_cell,'un',0));
|
|
|
+else
|
|
|
+ OutputMat = cell2mat(cellfun(@(x)cat(2,x',nan(1,maxLength-length(x))),input_cell,'un',0));
|
|
|
+end
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+%---------- plotting modules ----------%
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
|
|
+
|
|
|
+function plot_sca(ca_data, cellnum, clusternum)
|
|
|
+
|
|
|
+% putting data together to make plotting easier
|
|
|
+para = ca_data.stimpara;
|
|
|
+sta = [ca_data.locXhigh.sta;ca_data.locX_low.sta;ca_data.locY_high.sta;ca_data.locY_low.sta]';
|
|
|
+statime = - [ca_data.locXhigh.time;ca_data.locX_low.time;ca_data.locY_high.time;ca_data.locY_low.time]'*1e3;
|
|
|
+nlx = [ca_data.locXhigh.nlx;ca_data.locX_low.nlx;ca_data.locY_high.nlx;ca_data.locY_low.nlx]';
|
|
|
+nly = [ca_data.locXhigh.nly;ca_data.locX_low.nly;ca_data.locY_high.nly;ca_data.locY_low.nly]';
|
|
|
+
|
|
|
+figure('pos',[600 50 900 950],'color','w');
|
|
|
+% get plots title and legend labels
|
|
|
+titr = {'Location X', 'Location Y'};
|
|
|
+if para.twofixed
|
|
|
+ l = {'high-low','low-low'};
|
|
|
+ s = 'Spatial contrast adaptation';
|
|
|
+else
|
|
|
+ l = {'high-low','low-high'};
|
|
|
+ s = 'Spatial contrast adaptation (alternating stimulus)';
|
|
|
+end
|
|
|
+
|
|
|
+% first plotting psths
|
|
|
+subplot(3,2,[1,2])
|
|
|
+bar(para.low_contduration+para.psthtime_high(1:end-1),ca_data.psth_high);
|
|
|
+hold on;
|
|
|
+bar(para.psthtime_low(1:end-1),ca_data.psth_low);
|
|
|
+pbaspect([3 1 1]); box off;
|
|
|
+xlabel('Time (sec)'); ylabel('Rate (Hz)');
|
|
|
+xlim([0 para.low_contduration+para.high_contduration]);
|
|
|
+xticks([0 para.low_contduration para.low_contduration+para.high_contduration]);
|
|
|
+yAx = get(gca,'YLim'); yAx = ceil(yAx(2)/2)*2; yticks([0 yAx/2 yAx]);
|
|
|
+legend(l,'Location','northwest'); legend boxoff;
|
|
|
+
|
|
|
+% plotting STAs and nonlinearities
|
|
|
+for ii = 1:2
|
|
|
+
|
|
|
+ subplot(3,2,ii+2)
|
|
|
+ plot(statime(:,2*(ii-1)+1),sta(:,2*(ii-1)+1),statime(:,2*ii),sta(:,2*ii),'LineWidth',2);
|
|
|
+ pbaspect([1.2 1 1]); box off;
|
|
|
+ axis([0 600 -0.5 0.5]);
|
|
|
+ title(['STA ',titr{ii}]);
|
|
|
+ xlabel('Delay (sec)'); ylabel('Filter (a.u.)');
|
|
|
+ if ii==2, legend(l,'Location','southeast'); legend boxoff; end
|
|
|
+
|
|
|
+ subplot(3,2,ii+4)
|
|
|
+ plot(nlx(:,2*(ii-1)+1),nly(:,2*(ii-1)+1),nlx(:,2*ii),nly(:,2*ii),'LineWidth',2);
|
|
|
+ pbaspect([1.2 1 1]); box off; ax = gca;
|
|
|
+ yAx = ceil(ax.YLim(2)/2)*2;
|
|
|
+ axis([-3.1 3.1 0 yAx]);
|
|
|
+ title(['Nonlinearity ',titr{ii}]);
|
|
|
+ xlabel('Input'); ylabel('\Delta Rate (Hz)');
|
|
|
+ %ax.XAxisLocation = 'origin'; ax.YAxisLocation = 'origin';
|
|
|
+ yticks(0:yAx/2:yAx);
|
|
|
+ xticks(-3:1:3);
|
|
|
+
|
|
|
+end
|
|
|
+% mega title on top
|
|
|
+annotation('textbox', [0.1, 0.95, 0.8, 0], 'string', [s, ' properties of cell ',...
|
|
|
+ num2str(cellnum),', cluster ',num2str(clusternum)],'FontSize',16,'EdgeColor','none',...
|
|
|
+ 'HorizontalAlignment','center','VerticalAlignment','middle');
|
|
|
+
|
|
|
+end
|
|
|
+%--------------------------------------------------------------------------------------------------%
|
