% Anton Sumser 1st August 2022

load('ArticleData.mat','ROI_params','Trial_dFF','Trial_directionIdx','Directions','Trial_time','Trial_position','ROI_params_key','cell_count');
%% calculate direction selectivity

DSI=zeros(size(Trial_dFF,1),1); %initialize
DSIp=zeros(size(Trial_dFF,1),1); %initialize
prefDir=zeros(size(DSI,1),1);%initialize
norm_act=zeros(size(DSI,1),numel(Directions));
N_shuffles=1000;
for n=1:size(Trial_dFF,1)%run through ROIs
    Trial_Mean_act=mean(Trial_dFF{n}(:,Trial_position>0),2);
    [DSI(n),prefDir(n),norm_act(n,:)]=calculate_vector_dsi(Trial_Mean_act,Trial_directionIdx{n},Directions);
    DSIshuf=nan(1,N_shuffles);
    for s=1:N_shuffles
        DSIshuf(s)=calculate_vector_dsi(Trial_Mean_act,Trial_directionIdx{n}(randperm(size(Trial_directionIdx{n},1)),1),Directions);%calculate DSI with shuffled direction labels
    end
    DSIp(n)=sum(DSIshuf>repmat(DSI(n),1,N_shuffles),2)/N_shuffles;

end

%% Figure 6c

Varnames=cell_count.Properties.VariableNames;
figure
for a=1:3
plot(cell_count{:,'day of rec'},cell_count{:,Varnames{a+1}},'.');hold on;
end
cell_count_a=table2array(cell_count);
cell_count_a=cell_count_a(:,2:end);
plot(cell_count{:,'day of rec'},mean(cell_count_a,2,'omitnan'),'k');hold on

box off
xlim([2 17])
legend('animal 1','animal 2','animal 3','avg','location','northwest')
legend('boxoff')

%% Figure 6d

examples=[109 117 112 96 352 249   351  238];
tdtomato_examples=[-1 -1 -1 1 -1 -1 -1 1];
day_examples=[9 9 9 9 16 16 16 16];
directionsR=deg2rad(Directions);
mx=max(cat(1,Trial_dFF{examples}),[],'all');
figure('Position',[58          79        1182         899]);
tiledlayout(numel(examples),9);
for nx=1:numel(examples)
    n=examples(nx);
    
    for d=1:8
        nexttile;
        plot(Trial_time,Trial_dFF{n}(Trial_directionIdx{n}==d,:),'color',[.75 .75 .75],'LineWidth',.5);hold on;
        plot(Trial_time,mean(Trial_dFF{n}(Trial_directionIdx{n}==d,:),1),'k','LineWidth',1.5);hold off;box off
        ylim([0 6])
    end
    nexttile;
    pp=polarplot([directionsR;directionsR(1)],norm_act(n,[1:8 1]));
    rlim([0 .55])
    set(gca,'RTick',.25,'ThetaTick',[0 90 180 270])
end
%


%% Figure 6e
figure
[~,b]=histcounts(DSI,25);%bins for all ROIs
selectROIs=ROI_params(:,4)<.6; %tdTomato negative
DSI_select=DSI(selectROIs);
DSIp_select=DSIp(selectROIs);
[Ha]=histcounts(DSI_select,b);
[Hp05]=histcounts(DSI_select(DSIp_select<.05),b);
[Hp01]=histcounts(DSI_select(DSIp_select<.01),b);
bp=b(1:end-1)+mean(diff(b))/2;
bar(bp,Ha,1,'FaceAlpha',.5,'EdgeColor','none');hold on
bar(bp,Hp05,1,'FaceAlpha',.5,'EdgeColor','none');hold on
bar(bp,Hp01,1,'FaceAlpha',.5,'EdgeColor','none');hold on
xlabel('DSI');ylabel('#neurons')
legend('all p','p<.05','p<.01');legend('boxoff')
box off;
title('DSI histogram','tdTomato negative')

%% Figure 6f

figure
selectROIs=ROI_params(:,4)>=.6; %tdTomato positive
DSI_select=DSI(selectROIs);
DSIp_select=DSIp(selectROIs);
[Ha]=histcounts(DSI_select,b);
[Hp05]=histcounts(DSI_select(DSIp_select<.05),b);
[Hp01]=histcounts(DSI_select(DSIp_select<.01),b);

bar(bp,Ha,1,'FaceAlpha',.5,'EdgeColor','none');hold on
bar(bp,Hp05,1,'FaceAlpha',.5,'EdgeColor','none');hold on
bar(bp,Hp01,1,'FaceAlpha',.5,'EdgeColor','none');hold on
xlabel('DSI');ylabel('#neurons')
legend('all p','p<.05','p<.01');legend('boxoff')
box off;
title('DSI histogram','tdTomato positive')

%% Figure 6g
figure
polarscatter(prefDir(DSIp>=.05),DSI(DSIp>=.05),25,[0 0.4470 0.7410],'.');hold on
polarscatter(prefDir(DSIp<.05 & DSIp>=.01),DSI(DSIp<.05 & DSIp>=.01),25,[0.8500 0.3250 0.0980],'.');
polarscatter(prefDir(DSIp<.01),DSI(DSIp<.01),25,[0.9290 0.6940 0.1250],'.');hold off
legend('p>.05','p<.05','p<.01');legend('boxoff')
box off;
title('preferred directions')


%% Figure 6h
recording_day=ROI_params(:,6);
prop_sigDSI_day=nan(1,max(recording_day));
mean_DSI_day=nan(1,max(recording_day));
day_label=nan(1,max(recording_day));
for dy=1:max(recording_day)
    select_ROIs_on_day=find(recording_day==dy);
    prop_sigDSI_day(dy)=mean(DSIp(select_ROIs_on_day)<.05);
    mean_DSI_day(dy)=mean(DSI(select_ROIs_on_day));
    day_label(dy)=dy;
end
no_data=isnan(prop_sigDSI_day);
prop_sigDSI_day(no_data)=[];
day_label(no_data)=[];
mean_DSI_day(no_data)=[];

figure
plot(day_label,mean_DSI_day,'k');hold on
ylabel('mean DSI')
yyaxis right
plot(day_label,prop_sigDSI_day);
ylabel('proportion p(DSI)<.05')
box off
xlabel('days post CVS injection')


%% EDFigure 9c
figure
off=[-.25 0 .25];%fixed offset for plotting
for a=1:3
plot(ROI_params(ROI_params(:,5)==a,6)+off(a),ROI_params(ROI_params(:,5)==a,7),'.');hold on;%plor SNR over days per animal
end
legend('animal 1','animal 2','animal 3','Location','EastOutside');legend('boxoff')
xlabel('days post CVS injection')
ylabel('Neuron SNR')
box off




%% %%%%%%%%

%% helper functions
function [DSI,prefDir,norm_act]=calculate_vector_dsi(trialMeanData,directionIndices,directions)
%calculates direction selectivity index and preferred direction based on
%the "vector sum method":
%Mazurek, M., Kager, M. & Van Hooser, S. D. Robust quantification of orientation selectivity and direction selectivity. Frontiers in Neural Circuits 8, 92 (2014).

%INPUTS:
%trialMeanData: Mean calcium activity of a single neuron per trial (Trials x 1)
%directionIndices: index of direction per trial (Trials x 1) (can be shuffled)
%directions: Direction of movement directionIndices refer to in degrees (1 x max(directionIndices))

%OUTPUTS:
%DSI: direction selectivity index
%prefDir: preferred direction in rad

Direction_Average=nan(numel(directions),1);%initialize
for d=1:numel(directions)
    Direction_Average(d)=mean(trialMeanData(directionIndices==d)); %average trials with same direction
end
DSI=abs(sum(Direction_Average.*exp(1i*deg2rad(directions)),1))./sum(Direction_Average,1); % calculate direction selectivity
if nargout>1
    prefDir=angle(sum(Direction_Average.*exp(1i*deg2rad(directions)),1)); %calculate preferred directions
    norm_act=Direction_Average./sum(Direction_Average);
end
end