%decodes trial identity from spiking across bins for each individual neuron clear all load('RAWvppsALL.mat') load('RAWvpppsALL.mat') load('RAWvppasALL.mat') load('RAWvpppasALL.mat') if ~exist('RPPSall'),load('RPPSall.mat'); end if ~exist('RPSall'),load('RPSall.mat'); end if ~exist('RvppasALL'),load('RvppasALL.mat'); end if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end for PavInst=1:4; %DS is 1, Pav is 2 if PavInst==1 RAW=RAWvpppsALL; selection=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7; else if PavInst==2 RAW=RAWvppsALL; selection=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7; else if PavInst==3 RAW=RAWvpppasALL; selection=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7; else if PavInst==4 RAW=RAWvppasALL; selection=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7; end end end end TotalNeurons = 0; NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration repetitions=20; %how many times to run the analysis folds = 5; %number of times cross-validated: Dura=[0 0.3]; %size of bin used for decoding: bins = 1; %number of bins binint = 0.3; %spacing of bins binstart = 0; %start time of first bin relative to event shuffs = 1; %number of shuffled models created for i=1:length(NumNeurons) PoolMdlAcc{i,1}=NaN(folds,bins); PoolMdlAccShuff{i,1}=NaN(folds,bins); end %total number of neurons in all sessions and events per session %for i=1:length(RAW) TotalNeurons=sum(selection); % Ev1perSess(i)=length(RAW(i).Erast{Ev1,1}); % Ev2perSess(i)=length(RAW(i).Erast{Ev2,1}); %end %figure out how many trials of each event can be used % Ev1s=min(Ev1perSess); % Ev2s=min(Ev2perSess); %Need to keep number of trials used in VP and NAc constant, so have to pick %minimum number across all sessions in both regions. Right now with %existing sessions it's 20 Ev1s and 20 Ev2s Ev1s=27; Ev2s=27; for v = 1:length(NumNeurons) for u=1:repetitions %pick which neurons to use SetupSel=cat(1,ones(NumNeurons(v),1),zeros(TotalNeurons-NumNeurons(v),1)); Sel=(SetupSel(randperm(length(SetupSel)))==1); %setup the event identity matrix for decoding DecodeRs(1:Ev1s,1)=1; DecodeRs((Ev1s+1):(Ev1s+Ev2s),1)=2; for l=1:bins DecodeSpikes=NaN(1,1); AllSpikes=NaN(1,1); LowHz=zeros(1,NumNeurons(v)); NN = 1; AllNN=1; for i=1:length(RAW) %loops through sessions for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session if selection(AllNN)==1 Ev1Spikes=NaN(1,1); Ev2Spikes=NaN(1,1); %events being compared if (PavInst==1 | PavInst==2) Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2)); Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2)); else Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2)); Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2)); end %get number of spikes for this neuron in this bin for all %Ev1 trials [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials) for m=1:length(PSR1) Ev1Spikes(m,1)=sum(PSR1{1,m}>(binstart)); end %pick which trials get used for decoding SetupTrials=cat(1,ones(Ev1s,1),zeros(length(PSR1)-Ev1s,1)); Trials=SetupTrials==1; %(SetupTrials(randperm(length(SetupTrials)))==1); %randomize trials (or not) %put spikes from those trials in a matrix AllSpikes(Trials,NN)=Ev1Spikes(Trials,1); %get all the spikes from reward 1p2 trials [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials) for n=1:length(PSR2) Ev2Spikes(n,1)=sum(PSR2{1,n}>(binstart)); end %pick which trials get used for decoding SetupTrials2=cat(1,ones(Ev2s,1),zeros(length(PSR2)-Ev2s,1)); Trials2=SetupTrials2==1;%(SetupTrials2(randperm(length(SetupTrials2)))==1); %randomize trials (or not) %put spikes from those trials in a matrix AllSpikes((Ev1s+1):(Ev1s+Ev2s),NN)=Ev2Spikes(Trials2,1); NN=NN+1; %neuron counter AllNN=AllNN+1; else AllNN=AllNN+1; end end end %Setup spikes matrix for decoding DecodeSpikes=AllSpikes(:,Sel); %find neurons with too few spikes for t=1:NumNeurons(v) if sum(DecodeSpikes(1:Ev1s,t))<3 || sum(DecodeSpikes((1+Ev1s):(Ev1s+Ev2s),t))<3 LowHz(1,t)=1; end end %remove neurons with too few spikes % if sum(LowHz)>0 % Sel2=LowHz<1; % DecodeSpikes=DecodeSpikes(:,Sel2); % end %set up validation result matrices CVacc = NaN(folds,1); CVaccSh = NaN(folds,1); %normal model Partitions = cvpartition(DecodeRs,'KFold',folds); for r = 1:folds SVMModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)),'discrimType','diaglinear'); prediction = predict(SVMModel,DecodeSpikes(Partitions.test(r),:)); actual = DecodeRs(Partitions.test(r)); correct = prediction - actual; CVacc(r) = sum(correct==0) / length(correct); end %shuffled model for q=1:shuffs DecodeRsSh=DecodeRs(randperm(length(DecodeRs))); PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds); for s = 1:folds SVMModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'discrimType','diaglinear'); predictionSh = predict(SVMModelSh,DecodeSpikes(PartitionsSh.test(s),:)); actualSh = DecodeRsSh(PartitionsSh.test(s)); correctSh = predictionSh - actualSh; CVaccSh(s) = sum(correctSh==0) / length(correctSh); end AccShuff(q,1) = nanmean(CVaccSh); end PoolMdlAcc{v,1}(u,l) = nanmean(CVacc); PoolMdlAccShuff{v,1}(u,l) = nanmean(AccShuff); fprintf(['Condition #' num2str(v) ', Rep #' num2str(u) ', Bin #' num2str(l) '\n']); end %bins end end if PavInst==1 PoolMdlAccDSSucrose=PoolMdlAcc; PoolMdlAccDSSucroseShuff=PoolMdlAccShuff; save('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose') save('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff') else if PavInst==2 PoolMdlAccPavSucrose=PoolMdlAcc; PoolMdlAccPavSucroseShuff=PoolMdlAccShuff; save('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose') save('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff') else if PavInst==3 PoolMdlAccDSAlcohol=PoolMdlAcc; PoolMdlAccDSAlcoholShuff=PoolMdlAccShuff; save('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol') save('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff') else if PavInst==4 PoolMdlAccPavAlcohol=PoolMdlAcc; PoolMdlAccPavAlcoholShuff=PoolMdlAccShuff; save('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol') save('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff') end end end end end %% % % for PavInst=1:4; % if PavInst==1 % PoolMdlAccDSSucroseQuadratic=PoolMdlAccDSSucrose; % PoolMdlAccDSSucroseShuffQuadratic=PoolMdlAccDSSucroseShuff; % save('PoolMdlAccDSSucroseQuadratic.mat','PoolMdlAccDSSucroseQuadratic') % save('PoolMdlAccDSSucroseShuffQuadratic.mat','PoolMdlAccDSSucroseShuffQuadratic') % else if PavInst==2 % PoolMdlAccPavSucroseQuadratic=PoolMdlAccPavSucrose; % PoolMdlAccPavSucroseShuffQuadratic=PoolMdlAccPavSucroseShuff; % save('PoolMdlAccPavSucroseQuadratic.mat','PoolMdlAccPavSucroseQuadratic') % save('PoolMdlAccPavSucroseShuffQuadratic.mat','PoolMdlAccPavSucroseShuffQuadratic') % else if PavInst==3 % PoolMdlAccDSAlcoholQuadratic=PoolMdlAccDSAlcohol; % PoolMdlAccDSAlcoholShuffQuadratic=PoolMdlAccDSAlcoholShuff; % save('PoolMdlAccDSAlcoholQuadratic.mat','PoolMdlAccDSAlcoholQuadratic') % save('PoolMdlAccDSAlcoholShuffQuadratic.mat','PoolMdlAccDSAlcoholQuadratic') % else if PavInst==4 % PoolMdlAccPavAlcoholQuadratic=PoolMdlAccPavAlcohol; % PoolMdlAccPavAlcoholShuffQuadratic=PoolMdlAccPavAlcoholShuff; % save('PoolMdlAccPavAlcoholQuadratic.mat','PoolMdlAccPavAlcoholQuadratic') % save('PoolMdlAccPavAlcoholShuffQuadratic.mat','PoolMdlAccPavAlcoholShuffQuadratic') % % end % end % end % end % end