|
@@ -1,142 +0,0 @@
|
|
|
-function [pred_output]=BC_pred(STA_output)
|
|
|
-%[pred_output]=BC_pred(STA_output)
|
|
|
-%this function shows the prinicpal calculation for the the prediction to
|
|
|
-%the white noise stimulus with the linear-nonlinear model (LN-model)for
|
|
|
-%continous voltage signals of bipolar cells (BC).
|
|
|
-%Inputs:
|
|
|
-% STA_output = output of the BC_STA function; containing the
|
|
|
-% STA and stimulus signal
|
|
|
-%
|
|
|
-%Outpus: pred_output = structure containing the output of the
|
|
|
-% prediction
|
|
|
-%
|
|
|
-%
|
|
|
-% 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
|
|
|
-
|
|
|
-%-------------------------------------------------------------------------------
|
|
|
-%% 1. recompute STA and nonlinearity without the stimulus signal that we want to predict
|
|
|
-%%Important we want to predict a stimulus sequence that was not used to
|
|
|
-%%build the model (filter and nonlinearity). Therefore we have to redo the
|
|
|
-%%filter(STA) and nonlinearity without the stimulus-sequence that we want
|
|
|
-%to predict the response(see Methods of Schreyer&Gollisch,2021).
|
|
|
-%If we recorded Frozen frames, we predict the response to the frozen
|
|
|
-%frames from the model (filter and nonlinearity) we build from the running
|
|
|
-%frames.
|
|
|
-
|
|
|
-%% 1.1 STA
|
|
|
-%E.g. to make it easy: we want to predict the response to the last 300frames of the white
|
|
|
-%noise stimulus
|
|
|
-%take out the stimulus sequence you want to predict
|
|
|
-stimPred=STA_output.stimulus2D_zoom(:,end-299:end); %the last 300frames
|
|
|
-%voltage trace to predict
|
|
|
-voltPred=STA_output.membranPotAVG(end-299:end); %the response to the last 300 frames
|
|
|
-%build the STA without the last 300 frames
|
|
|
-[STA_norm]=buildSTA ...
|
|
|
- (STA_output.stimulus2D_zoom(:,1:end-300),...
|
|
|
- STA_output.membranPotAVG(1:end-300),STA_output.STA_fnbr);
|
|
|
-STA3D_norm=reshape(STA_norm,size(STA_output.rangeY_zoom,2),size(STA_output.rangeX_zoom,2), ...
|
|
|
- size(STA_norm,2));% Y-zoom checkers (monitor axis), X-zoom checkers (monitor axis) and Z (STA_fnbr)
|
|
|
-
|
|
|
-%% 1.2 Nonlinearity
|
|
|
-%build the NL from the data without the last 300 frames
|
|
|
-generator_signal=filter2(STA_norm,STA_output.stimulus2D_zoom(:,1:end-300),'valid');
|
|
|
-%get the corresponding membrane potential responses
|
|
|
-y_axisNL_unsorted=STA_output.membranPotAVG(STA_output.STA_fnbr:end-300); %the first value is at the length of the filter
|
|
|
-
|
|
|
-%sort the x axis
|
|
|
-[NL_xAxis,indx]=sort(generator_signal);
|
|
|
-%sort the corresponding y axis
|
|
|
-NL_yAxis=y_axisNL_unsorted(indx);
|
|
|
-%do a binning with percentile to have same amount of data in each bin
|
|
|
-%make 40 bins
|
|
|
-nbrbins=40;
|
|
|
-bins= doBin(NL_xAxis,NL_yAxis,nbrbins);
|
|
|
-
|
|
|
-%-------------------------------------------------------------------------------
|
|
|
-%% 2. compute the prediction
|
|
|
-%% 2.1 convolution
|
|
|
-%convolve the new stimulus sequence (not used for prediciton) with the
|
|
|
-%normalized STA
|
|
|
-stim2pred=filter2(STA_norm,stimPred,'valid');
|
|
|
-
|
|
|
-%% 2.2 compute the prediction through linear interpolation
|
|
|
-NL_x=bins.NL_xbin;
|
|
|
-NL_y=bins.NL_ybin;
|
|
|
-%example of how to make the prediction through linear interpolation
|
|
|
-predMemP_y=nan(1,size(stim2pred,2));
|
|
|
-for kk=1:size(stim2pred,2)
|
|
|
- nbrTemp=stim2pred(kk); %first x-axis value to predict
|
|
|
- [xIndxS,xIndxL]= getIndex(NL_x,nbrTemp,1); %get those values that are on left and rigth
|
|
|
- xValueS=NL_x(xIndxS);
|
|
|
- xValueL=NL_x(xIndxL);
|
|
|
- %do the linear interpolation
|
|
|
- [polyCoeff] = polyfit([xValueS;xValueL],[NL_y(xIndxS);...
|
|
|
- NL_y(xIndxL)],1);
|
|
|
- predMemP_y(kk)=polyCoeff(1)*nbrTemp+polyCoeff(2);
|
|
|
-end
|
|
|
-
|
|
|
-%-------------------------------------------------------------------------------
|
|
|
-%% 3. evaluate the prediction performance
|
|
|
-%pearson correlation
|
|
|
-corrcoeff_temp=corrcoef(voltPred(STA_output.STA_fnbr:end),predMemP_y); %the first value is at the length of the filter
|
|
|
-pearCorrSqrt=corrcoeff_temp(2)^2;
|
|
|
-
|
|
|
-%-------------------------------------------------------------------------------
|
|
|
-%% 4. plot the output
|
|
|
-figure;
|
|
|
-%STA
|
|
|
-subplot(2,2,1)
|
|
|
-[~,indx_Zoom]=max(abs(STA3D_norm(:))); %max intensity value of STA
|
|
|
-[~,~,posZ]=ind2sub(size(STA3D_norm),indx_Zoom); %get the frame position
|
|
|
-imagesc(STA3D_norm(:,:,posZ))
|
|
|
-colormap('gray'); %how the sitmulus was shown
|
|
|
-axis equal; axis tight;
|
|
|
-title('STA best frame')
|
|
|
-%NL
|
|
|
-subplot(2,2,2)
|
|
|
-plot(NL_xAxis,NL_yAxis,'.','color','k') %plot non-binned signal
|
|
|
-hold on;
|
|
|
-plot(bins.NL_xbin,bins.NL_ybin,'o','markerfacecolor','r','markeredgecolor','r') %plot binned signal
|
|
|
-title('NL')
|
|
|
-axis tight
|
|
|
-xlabel('generator signal')
|
|
|
-ylabel('voltage (mV)')
|
|
|
-%prediction
|
|
|
-subplot(2,2,3:4)
|
|
|
-plot(voltPred(STA_output.STA_fnbr:end),'k') %original membrane potential
|
|
|
-hold on;
|
|
|
-plot(predMemP_y,'r')%predicted membrane potential
|
|
|
-title(['pCorr:' mat2str(round(pearCorrSqrt,2)) ' Filename: ' STA_output.filename])
|
|
|
-xlabel('stimulus frames')
|
|
|
-ylabel('voltage (ms)')
|
|
|
-legend('original response','predicted response')
|
|
|
-
|
|
|
-%-------------------------------------------------------------------------------
|
|
|
-%% 5. add variables into the output
|
|
|
-pred_output.pearCorrSqrt=pearCorrSqrt;
|
|
|
-
|
|
|
-end
|
|
|
-
|
|
|
-function [IndS,IndL]= getIndex(x,valConv,Points)
|
|
|
-IndS=find(x<valConv,Points,'last'); %find the last 5 closest smaller values
|
|
|
-IndL=find(x>=valConv,Points);
|
|
|
-if isempty(IndS)|| length(IndS)<Points %if you are on the left border
|
|
|
- if isempty(IndS) %take value from 1:10
|
|
|
- IndL=find(x>=valConv,Points*2);
|
|
|
- IndS=IndL(1);
|
|
|
- IndL=IndL(2:end);
|
|
|
- else
|
|
|
- IndL=find(x>=valConv,Points*2-length(IndS));
|
|
|
- end
|
|
|
-elseif isempty(IndL)|| length(IndL)<Points %if you are on the right border
|
|
|
- if isempty(IndL) %take value from 1:10
|
|
|
- IndS=find(x<valConv,Points*2,'last');
|
|
|
- IndL=IndS(end);
|
|
|
- IndS=IndS(1:end-1);
|
|
|
- else
|
|
|
- IndS=find(x<valConv,Points*2-length(IndL));
|
|
|
- end
|
|
|
-end
|
|
|
-end
|