clear;
close all; 
% paired analysis per unit tested in different protocols
load('pair_data.mat')

class_=all_di(:,2)==0;  % unit classification  
paired_=logical([1;1;1;1;1;1;0;0;0;0;0;0;0;0;0;1;1;0;1;1;1;1;1;1;0;0;1;1]);
aw_=all_di(:,3)==0;  % awake data 
an_=all_di(:,3)==1;  % aneshetized data 
particular=logical(aw_);
an_units=all_di(particular,1)+1;
particular2=ismember(all_di(:,1),an_units);
n=0;  % 0: for 60ms delay and 3: for 416ms delay
c=0;  % 0: for 60ms delay and 4: for 416ms delay

% respective discriminability index
di_iso=all_di(particular,4+n);
di_nav=all_di(particular,5+n);   
di_com=all_di(particular,6+n); 

di_iso_comp=all_di(particular2,4+n);
di_nav_comp=all_di(particular2,5+n);
di_com_comp=all_di(particular2,6+n);

y1=[di_nav,di_iso,di_iso_comp,di_nav_comp];
y2=[di_com,di_iso,di_iso_comp,di_com_comp];

% respective context effect index
ce_nav_exp=all_di(particular,10+c);
ce_nav_unexp=all_di(particular,11+c);
ce_com_exp=all_di(particular,12+c);
ce_com_unexp=all_di(particular,13+c);

ce_nav_exp_comp=all_di(particular2,10+c);
ce_nav_unexp_comp=all_di(particular2,11+c);
ce_com_exp_comp=all_di(particular2,12+c);
ce_com_unexp_comp=all_di(particular2,13+c);

y3=[ce_nav_exp,ce_nav_exp_comp,ce_nav_unexp,ce_nav_unexp_comp];
y4=[ce_com_exp,ce_com_exp_comp,ce_com_unexp,ce_com_unexp_comp];

% respective deviance detection index
dd_nav=ce_nav_unexp-ce_nav_exp;
dd_com=ce_com_unexp-ce_com_exp;

dd_nav_comp=ce_nav_unexp_comp-ce_nav_exp_comp;
dd_com_comp=ce_com_unexp_comp-ce_com_exp_comp;

y5=[dd_nav,dd_nav_comp,dd_com,dd_com_comp];

% respective spikecounts in isolation
sp_nav=all_di(particular,18);
sp_com=all_di(particular,19);

sp_nav_comp=all_di(particular2,18);
sp_com_comp=all_di(particular2,19);

y6=[sp_nav,sp_nav_comp,sp_com,sp_com_comp];

% compare paired samples discriminability index 
pvA1=signrank(di_iso,di_nav);
pvA2=signrank(di_iso_comp,di_nav_comp);
pvA3=signrank(di_iso,di_com);
pvA4=signrank(di_iso_comp,di_com_comp);
pvA5=signrank(di_nav,di_nav_comp);
pvA6=signrank(di_com,di_com_comp);

% compare paired samples context effect index 
pvB1=signrank(ce_nav_exp,ce_nav_exp_comp);
pvB2=signrank(ce_nav_unexp,ce_nav_unexp_comp);
pvB3=signrank(ce_nav_exp,ce_nav_unexp);
pvB4=signrank(ce_nav_exp_comp,ce_nav_unexp_comp);
pvB5=signrank(ce_com_exp,ce_com_exp_comp);
pvB6=signrank(ce_com_unexp,ce_com_unexp_comp);
pvB7=signrank(ce_com_exp,ce_com_unexp);
pvB8=signrank(ce_com_exp_comp,ce_com_unexp_comp);

% compare paired samples discriminability index 
pvC1=signrank(dd_nav,dd_nav_comp);
pvC2=signrank(dd_com,dd_com_comp);

% compare paired samples spikecounts
pvD1=signrank(sp_nav,sp_nav_comp);
pvD2=signrank(sp_com,sp_com_comp);


%% plot d.i. in isolation and after two different context  
red =[0.8500,0.3250, 0.0980];
blue=[0,0.4470,0.7410];
red_light=[1,0.5,0.5].*red;
blue_light=[1,0.5,0.5].*blue;
black=[0.2 0.2 0.2];

lines_color=[0.7 0.7 0.7];  % light gray

bar_colors=[blue_light; black; black; blue];
x = 1:4;

figure(1); set(gcf,'Position',[400 200 200 250])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y1(:,i))), y1(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end

% add lines per unit
hold on
for i=1:numel(di_iso)
    l=plot([1,2], [di_nav(i), di_iso(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [di_iso_comp(i), di_nav_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylabel('discriminability index')
xticklabels({'bef-nav','bef-sil','aft-sil','aft-nav'});


set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=0.9;
line2=1.1;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot([xt(1) xt(4)], [1 1]*line2, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line2,line2-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line2,line2-0.05],'-k', 'LineWidth',sl);


ylim([-1 1.3])

% plot level of significance
if pvA1 <= 0.05 && pvA1 > 0.01
    text(mean(xt(1:2))-0.025, line1+0.05, '*')
elseif pvA1 <=0.01 && pvA1 > 0.001
    text(mean(xt(1:2))-0.05, line1+0.05, '**')
elseif pvA1 <= 0.001
    text(mean(xt(1:2))-0.15, line1+0.05, '***')
else 
    text(mean(xt(1:2))-0.05, line1+0.075, 'ns','FontSize',8)
end

if pvA2 <= 0.05 && pvA2 > 0.01
    text(mean(xt(3:4))-0.025, line12+0.05, '*')
elseif pvA2 <=0.01 && pvA2 > 0.001
    text(mean(xt(3:4))-0.05, line1+0.05, '**')
elseif pvA2 <= 0.001
    text(mean(xt(3:4))-0.15, line1+0.05, '***')
else 
    text(mean(xt(3:4))-0.05, line1+0.075, 'ns','FontSize',8)
end

if pvA5 <= 0.05 && pvA5 > 0.01
    text(mean(xt(1:4))-0.025, line2+0.05, '*')
elseif pvA5 <=0.01 && pvA5 > 0.001
    text(mean(xt(1:4))-0.05, line21+0.05, '**')
elseif pvA5 <= 0.001
    text(mean(xt(1:4))-0.15, line2+0.05, '***')
else 
    text(mean(xt(1:4))-0.05, line2+0.075, 'ns','FontSize',8)
end


bar_colors=[red_light; black; black; red];
lines_color=[0.7 0.7 0.7];  % light gray
x = 1:4;

figure(2); set(gcf,'Position',[600 200 200 250])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y2(:,i))), y2(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end

% add lines per unit
hold on
for i=1:numel(di_iso)
    l=plot([1,2], [di_com(i), di_iso(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [di_iso(i), di_com_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylabel('discriminability index')
xticklabels({'bef-com','bef-sil','aft-sil','aft-com'});
set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=0.9;
line2=1.1;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot([xt(1) xt(4)], [1 1]*line2, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line2,line2-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line2,line2-0.05],'-k', 'LineWidth',sl);


ylim([-1 1.3])

% plot level of significance
if pvA3 <= 0.05 && pvA3 > 0.01
    text(mean(xt(1:2))-0.025, line1+0.05, '*')
elseif pvA3 <=0.01 && pvA3 > 0.001
    text(mean(xt(1:2))-0.05, line1+0.05, '**')
elseif pvA3 <= 0.001
    text(mean(xt(1:2))-0.15, line1+0.05, '***')
else 
    text(mean(xt(1:2))-0.05, line1+0.075, 'ns','FontSize',8)
end

if pvA4 <= 0.05 && pvA4 > 0.01
    text(mean(xt(3:4))-0.025, line1+0.05, '*')
elseif pvA4 <=0.01 && pvA4 > 0.001
    text(mean(xt(3:4))-0.05, line1+0.05, '**')
elseif pvA4 <= 0.001
    text(mean(xt(3:4))-0.15, line1+0.05, '***')
else 
    text(mean(xt(3:4))-0.05, line1+0.075, 'ns','FontSize',8)
end

if pvA6 <= 0.05 && pvA6 > 0.01
    text(mean(xt(1:4))-0.025, line2+0.05, '*')
elseif pvA6 <=0.01 && pvA6 > 0.001
    text(mean(xt(1:4))-0.05, line21+0.05, '**')
elseif pvA6 <= 0.001
    text(mean(xt(1:4))-0.15, line2+0.05, '***')
else 
    text(mean(xt(1:4))-0.05, line2+0.075, 'ns','FontSize',8)
end

%% plot c.e. for expected and unexpected probes after two different context  
bar_colors=[blue;red;blue_light;red_light];

x = 1:4;

figure(3); set(gcf,'Position',[400 200 200 200])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y3(:,i))), y3(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(ce_nav_exp)
    l=plot([1,2], [ce_nav_exp(i), ce_nav_exp_comp(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [ce_nav_unexp(i), ce_nav_unexp_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylabel('echolocation effect')
xticklabels({'match-aw','match-an','mismatch-aw','mismatch-an'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=0.45;
line2=0.65;
line3=0.85;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.05],'-k', 'LineWidth',sl);

% plot([xt(1) xt(3)], [1 1]*line2, '-k','LineWidth',sl) 
% plot([1;1]*xt(1),[line2,line2-0.05],'-k', 'LineWidth',sl);
% plot([1;1]*xt(3),[line2,line2-0.05],'-k', 'LineWidth',sl);
% 
% plot([xt(2) xt(4)], [1 1]*line3, '-k','LineWidth',sl) 
% plot([1;1]*xt(2),[line3,line3-0.05],'-k', 'LineWidth',sl);
% plot([1;1]*xt(4),[line3,line3-0.05],'-k', 'LineWidth',sl);

ylim([-1.05 0.7])

% % plot level of significance
% if pvB1 <= 0.05 && pvB1 > 0.01
%     text(mean(xt(1:2))-0.025, line1+0.05, '*')
% elseif pvB1 <=0.01 && pvB1 > 0.001
%     text(mean(xt(1:2))-0.05, line1+0.05, '**')
% elseif pvB1 <= 0.001
%     text(mean(xt(1:2))-0.15, line1+0.05, '***')
% else 
%     text(mean(xt(1:2))-0.05, line1+0.1, 'ns','FontSize',7)
% end
% 
% if pvB2 <= 0.05 && pvB2 > 0.01
%     text(mean(xt(3:4))-0.025, line1+0.05, '*')
% elseif pvB2 <=0.01 && pvB2 > 0.001
%     text(mean(xt(3:4))-0.05, line1+0.05, '**')
% elseif pvB2 <= 0.001
%     text(mean(xt(3:4))-0.15, line1+0.05, '***')
% else 
%     text(mean(xt(3:4))-0.05, line1+0.1, 'ns','FontSize',7)
% end
% 
% 
% if pvB3 <= 0.05 && pvB3 > 0.01
%     text(mean(xt(1:3))-0.025, line2+0.05, '*')
% elseif pvB3 <=0.01 && pvB3 > 0.001
%     text(mean(xt(1:3))-0.05, line2+0.05, '**')
% elseif pvB3 <= 0.001
%     text(mean(xt(1:3))-0.15, line2+0.05, '***')
% else 
%     text(mean(xt(1:3))-0.05, line2+0.1, 'ns','FontSize',7)
% end
% 
% 
% if pvB4 <= 0.05 && pvB4 > 0.01
%     text(mean(xt(2:4))-0.025, line3+0.05, '*')
% elseif pvB4 <=0.01 && pvB4 > 0.001
%     text(mean(xt(2:4))-0.05, line3+0.05, '**')
% elseif pvB4 <= 0.001
%     text(mean(xt(2:4))-0.15, line3+0.05, '***')
% else 
%     text(mean(xt(2:4))-0.05, line3+0.1, 'ns','FontSize',7)
% end

exportgraphics(gcf,'E:\Users\User\Desktop\delay paper\anesthetized\figures_paper\sameunit_inj_nav.pdf',...
    'Resolution',300','ContentType','vector','BackgroundColor','none')


figure(4); set(gcf,'Position',[400 200 200 200])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y4(:,i))), y4(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(ce_nav_exp)
    l=plot([1,2], [ce_com_exp(i), ce_com_exp_comp(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [ce_com_unexp(i), ce_com_unexp_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylim([-1.05, 0.4])
ylabel('communication effect')
xticklabels({'match-aw','match-an','mismatch-aw','mismatch-an'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=0.45;
line2=0.65;
line3=0.85;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.05],'-k', 'LineWidth',sl);

% plot([xt(1) xt(3)], [1 1]*line2, '-k','LineWidth',sl) 
% plot([1;1]*xt(1),[line2,line2-0.05],'-k', 'LineWidth',sl);
% plot([1;1]*xt(3),[line2,line2-0.05],'-k', 'LineWidth',sl);
% 
% plot([xt(2) xt(4)], [1 1]*line3, '-k','LineWidth',sl) 
% plot([1;1]*xt(2),[line3,line3-0.05],'-k', 'LineWidth',sl);
% plot([1;1]*xt(4),[line3,line3-0.05],'-k', 'LineWidth',sl);
% 
ylim([-1.05 0.7])
% 
% % plot level of significance
% if pvB5 <= 0.05 && pvB5 > 0.01
%     text(mean(xt(1:2))-0.025, line1+0.05, '*')
% elseif pvB5 <=0.01 && pvB5 > 0.001
%     text(mean(xt(1:2))-0.05, line1+0.05, '**')
% elseif pvB5 <= 0.001
%     text(mean(xt(1:2))-0.15, line1+0.05, '***')
% else 
%     text(mean(xt(1:2))-0.05, line1+0.1, 'ns','FontSize',7)
% end
% 
% if pvB6 <= 0.05 && pvB6 > 0.01
%     text(mean(xt(3:4))-0.025, line1+0.05, '*')
% elseif pvB6 <=0.01 && pvB6 > 0.001
%     text(mean(xt(3:4))-0.05, line1+0.05, '**')
% elseif pvB6 <= 0.001
%     text(mean(xt(3:4))-0.15, line1+0.05, '***')
% else 
%     text(mean(xt(3:4))-0.05, line1+0.1, 'ns','FontSize',7)
% end
% 
% 
% if pvB7 <= 0.05 && pvB7 > 0.01
%     text(mean(xt(1:3))-0.025, line2+0.05, '*')
% elseif pvB7 <=0.01 && pvB7 > 0.001
%     text(mean(xt(1:3))-0.05, line2+0.05, '**')
% elseif pvB7 <= 0.001
%     text(mean(xt(1:3))-0.15, line2+0.05, '***')
% else 
%     text(mean(xt(1:3))-0.05, line2+0.1, 'ns','FontSize',7)
% end
% 
% 
% if pvB8 <= 0.05 && pvB8 > 0.01
%     text(mean(xt(2:4))-0.025, line3+0.05, '*')
% elseif pvB8 <=0.01 && pvB8 > 0.001
%     text(mean(xt(2:4))-0.05, line3+0.05, '**')
% elseif pvB8 <= 0.001
%     text(mean(xt(2:4))-0.15, line3+0.05, '***')
% else 
%     text(mean(xt(2:4))-0.05, line3+0.1, 'ns','FontSize',7)
% end

exportgraphics(gcf,'E:\Users\User\Desktop\delay paper\anesthetized\figures_paper\sameunit_inj_com.pdf',...
    'Resolution',300','ContentType','vector','BackgroundColor','none')


%% plot d.d. for different context after two different context  
bar_colors=[blue_light;blue;red_light;red];

x = 1:4;

figure(5); set(gcf,'Position',[400 200 250 250])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y5(:,i))), y5(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(dd_nav)
    l=plot([1,2], [dd_nav(i), dd_nav_comp(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [dd_com(i), dd_com_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylabel('deviance detection')
xticklabels({'before','after','before','after'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=0.75;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.05],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.05],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.05],'-k', 'LineWidth',sl);

ylim([-0.75 1])

% plot level of significance
if pvC1 <= 0.05 && pvC1 > 0.01
    text(mean(xt(1:2))-0.025, line1+0.05, '*')
elseif pvC1 <=0.01 && pvC1 > 0.001
    text(mean(xt(1:2))-0.05, line1+0.05, '**')
elseif pvC1 <= 0.001
    text(mean(xt(1:2))-0.15, line1+0.05, '***')
else 
    text(mean(xt(1:2))-0.05, line1+0.075, 'ns','FontSize',8)
end

if pvC2 <= 0.05 && pvC2 > 0.01
    text(mean(xt(3:4))-0.025, line1+0.05, '*')
elseif pvC2 <=0.01 && pvC2 > 0.001
    text(mean(xt(3:4))-0.05, line1+0.05, '**')
elseif pvC2 <= 0.001
    text(mean(xt(3:4))-0.15, line1+0.05, '***')
else 
    text(mean(xt(3:4))-0.05, line1+0.075, 'ns','FontSize',8)
end

%% plot spikecounts of each probe after two different context  
bar_colors=[blue_light;blue;red_light;red];

x = 1:4;

figure(6); set(gcf,'Position',[400 200 250 250])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y6(:,i))), y6(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(dd_nav)
    l=plot([1,2], [sp_nav(i), sp_nav_comp(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [sp_com(i), sp_com_comp(i)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 4.5])
ylabel('# spikes in 50 ms')
xticklabels({'before','after','before','after'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=15;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-0.5],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-0.5],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-0.5],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-0.5],'-k', 'LineWidth',sl);

ylim([0 17])

% plot level of significance
if pvD1 <= 0.05 && pvD1 > 0.01
    text(mean(xt(1:2))-0.025, line1+0.5, '*')
elseif pvD1 <=0.01 && pvD1 > 0.001
    text(mean(xt(1:2))-0.05, line1+0.5, '**')
elseif pvD1 <= 0.001
    text(mean(xt(1:2))-0.15, line1+0.5, '***')
else 
    text(mean(xt(1:2))-0.05, line1+0.75, 'ns','FontSize',8)
end

if pvD2 <= 0.05 && pvD2 > 0.01
    text(mean(xt(3:4))-0.025, line1+0.5, '*')
elseif pvD2 <=0.01 && pvD2 > 0.001
    text(mean(xt(3:4))-0.05, line1+0.5, '**')
elseif pvD2 <= 0.001
    text(mean(xt(3:4))-0.15, line1+0.5, '***')
else 
    text(mean(xt(3:4))-0.05, line1+0.75, 'ns','FontSize',8)
end