Schreyer_and_Gollisch_2021_salamander_retinal_bipolar_cell_recordings/spatiotemporalWhiteNoise/code spatiotemporal white noise/BC_pred.m

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 March 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).
%In the Manuscript we predict 200 sequences of 300 frames long that where
%randomly chosen (see Method section). Every time we rebuilt the STA and
%NL. The overall prediction performance is then the average R^2 over all
%200 sequences. Here, for simplicity I only show the prediction to one
%sequence. Note however, to get a good estimate of how well the LN-model can predict
%the response of the cell it is important to run it over multiple
%sequences.
%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
%to keep the code easy: I here show the principal idea, by predicting the 
%response to the last 300frames of the white noise stimulus. 
%(in the manuscript we randomly choose 200 sequences and redo this part 200 times)

%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'); %see also Point 7
%in the function BC_STA.m and BC_NL of how you do it when you have a best
%spatial and temporal component and reduce noice by cutting 3-sigma
%contour around the receptive field region.
%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