%% Figure 2

load('Figure 2.mat')

sessions_numb=10;

figure('units','normalized','outerposition',[0.2 0.4 0.38 0.38])
hold on

plot( 1:sessions_numb,100*cell2mat(CR_pecr_avg_ipsi(:,1:sessions_numb,2)),'color', 'b' )

standart_err=100*nanstd( cell2mat(CR_pecr_avg_ipsi(:,:,1) ) )./sqrt(sum(  isnan(cell2mat(CR_pecr_avg_ipsi(:,:,1)))==0  ));   
errorbar( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_ipsi(:,1:sessions_numb,1))),standart_err(1:sessions_numb),'color', 'k')
scatter( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_ipsi(:,1:sessions_numb,1))),70, 'k','filled')

standart_err=100*nanstd( cell2mat(CR_pecr_avg_ipsi(:,:,2)) )./sqrt(sum(  isnan(cell2mat(CR_pecr_avg_ipsi(:,:,2)))==0  ));   
errorbar( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_ipsi(:,1:sessions_numb,2))),standart_err(1:sessions_numb),'color', 'b')
scatter( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_ipsi(:,1:sessions_numb,2))),70, 'b', 'filled'  )

xlabel('Session number')
ylabel('CR probability')
set(gca,'YTick',0:20:100)
set(gca,'XTick',0:2:100)
set(gca,'LineWidth',1)
box off
axis([1 sessions_numb 0 100])


figure('units','normalized','outerposition',[0.2 0.4 0.38 0.38])
hold on

plot( 1:sessions_numb,cell2mat(CR_amp_avg_ipsi(:,1:sessions_numb,2)),'color', 'b')

standart_err=nanstd( cell2mat(CR_amp_avg_ipsi(:,:,1) ) )./sqrt(sum(  isnan(cell2mat(CR_amp_avg_ipsi(:,:,1)))==0  ));   
errorbar( 1:sessions_numb,nanmean(cell2mat(CR_amp_avg_ipsi(:,1:sessions_numb,1))),standart_err(1:sessions_numb),'color', 'k')
scatter( 1:sessions_numb,nanmean(cell2mat(CR_amp_avg_ipsi(:,1:sessions_numb,1))),70, 'k','filled')

standart_err=nanstd( cell2mat(CR_amp_avg_ipsi(:,:,2)) )./sqrt(sum(  isnan(cell2mat(CR_amp_avg_ipsi(:,:,2)))==0  ));   
errorbar( 1:sessions_numb,nanmean(cell2mat(CR_amp_avg_ipsi(:,1:sessions_numb,2))),standart_err(1:sessions_numb),'color', 'b')
scatter( 1:sessions_numb,nanmean(cell2mat(CR_amp_avg_ipsi(:,1:sessions_numb,2))),70, 'b', 'filled'  )

xlabel('Session number')
ylabel('Mean CR amplitude, mm')

set(gca,'YTick',0:2:6)
set(gca,'XTick',0:2:10)
set(gca,'LineWidth',1)
box off
axis([1 sessions_numb 0 6])

clear

%% Figure 3
load('Figure 3.mat')

sessions_numb=10;

figure('units','normalized','outerposition',[0.2 0.4 0.38 0.38])
hold on

plot( 1:sessions_numb,100*cell2mat(CR_pecr_avg_contra(:,1:sessions_numb,2)),'color', 'r'  )

standart_err=100*nanstd( cell2mat(CR_pecr_avg_contra(:,:,2) ) )./sqrt(sum(  isnan(cell2mat(CR_pecr_avg_contra(:,:,2) ))==0  ));   
errorbar( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_contra(:,1:sessions_numb,2))),standart_err(1:sessions_numb),'color', 'r')
scatter( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_contra(:,1:sessions_numb,2))),70, 'r','filled')

standart_err=100*nanstd( cell2mat(CR_pecr_avg_contra(:,:,1) ) )./sqrt(sum(  isnan(cell2mat(CR_pecr_avg_contra(:,:,1)))==0  ));   
errorbar( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_contra(:,1:sessions_numb,1))),standart_err(1:sessions_numb),'color', 'k')
scatter( 1:sessions_numb,100*nanmean(cell2mat(CR_pecr_avg_contra(:,1:sessions_numb,1))),70, 'k','filled')

xlabel('Session number')
ylabel('CR probability')
set(gca,'YTick',0:20:100)
set(gca,'XTick',0:2:10)
set(gca,'LineWidth',1)
box off
axis([1 10 0 100])

figure('units','normalized','outerposition',[0.2 0.4 0.38 0.38])
hold on
    
plot( 1:sessions_numb,cell2mat(CR_amp_contra(:,1:sessions_numb,2)),'color','r'  )

standart_err=nanstd( cell2mat(CR_amp_contra(:,:,2) ) )./sqrt(sum(  isnan(cell2mat(CR_amp_contra(:,:,2)))==0  ));   
errorbar( 1:sessions_numb,nanmean(cell2mat(CR_amp_contra(:,1:sessions_numb,2))),standart_err(1:sessions_numb),'color', 'r')
scatter( 1:sessions_numb,nanmean(cell2mat(CR_amp_contra(:,1:sessions_numb,2))),70, 'r','filled')

standart_err=nanstd( cell2mat(CR_amp_contra(:,:,1) ) )./sqrt(sum(  isnan(cell2mat(CR_amp_contra(:,:,1)))==0  ));   
errorbar( 1:sessions_numb,nanmean(cell2mat(CR_amp_contra(:,1:sessions_numb,1))),standart_err(1:sessions_numb),'color', 'k')
scatter( 1:sessions_numb,nanmean(cell2mat(CR_amp_contra(:,1:sessions_numb,1))),70, 'k','filled')


xlabel('Session number')
ylabel('Mean CR amplitude, mm')
set(gca,'XTick',0:2:10)
set(gca,'YTick',0:2:6)
set(gca,'LineWidth',1)
box off
axis([1 10 0 6])

clear

%% Figure 4

load('Figure 4.mat')

figure('units','normalized','outerposition',[0.2 0.4 0.5 0.5]) 
hold on
set(gca,'LineWidth',1)
box off
CR_amps_tmp=cell2mat(CR_amps_ipsi(:));

scatter(CR_amps_tmp(:,1),CR_amps_tmp(:,2),15,'k','filled')

plot([-1 10], [-1 10],'k')
plot([0 10], [0.3 0.3],'k')
plot([0.3 0.3],[0 10],'k')
axis equal
axis([-0. 7 -0. 7])
set(gca,'YTick',0:2:6)
set(gca,'XTick',0:2:6)
xlabel('First CR amplitude, mm')
ylabel('Second CR amplitude, mm')

% conditional probabiltity of 2nd CR  given 1st CR amp
amp_bins=[0 0.3 2.5:6.5];

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4]) 
hold on
set(gca,'LineWidth',1)
box off

first_CR_amp=CR_amps_tmp(:,1);
second_CR_amp=CR_amps_tmp(:,2);
amp_coord=[];
prob_2nd_CR_conf=[];
prob_2nd_CR=[];

for i=1:length(amp_bins)-1
    tr=find(first_CR_amp>amp_bins(i)&first_CR_amp<amp_bins(i+1));
    second_CR_amp_tr=second_CR_amp(tr);
    prob_2nd_CR(i)=sum(second_CR_amp_tr>0.3)/length(second_CR_amp_tr);
    amp_coord(i)=(amp_bins(i+1)+amp_bins(i))/2;
    prob_2nd_CR_conf(i,:)=bootci(1000,{@(x) (sum(x>0.3))/length(x),second_CR_amp_tr},'alpha',0.05);
end

scatter(amp_coord, 100*prob_2nd_CR,100, 'k' ,'filled')
plot(amp_coord, 100*prob_2nd_CR, 'k' )
errorbar( amp_coord, 100*prob_2nd_CR, 100*prob_2nd_CR-100*prob_2nd_CR_conf(:,1)',100*prob_2nd_CR_conf(:,2)'- 100*prob_2nd_CR, 'k','MarkerSize',1 ,'LineWidth',1)

axis([0 6 0 100])
xlabel('First CR amplitude, mm')
ylabel('Second CR probability')
set(gca,'YTick',0:20:100)
set(gca,'XTick',0:2:6)


% First and second CR amplitude during Ext. 1st Reinf. 2nd CR sessions 

block_number=6;
bunny_numb=4;
max_sessions_numb=2;

CR_amp_block=NaN(bunny_numb,max_sessions_numb*block_number,2);
CR_amp_block=mat2cell(CR_amp_block,ones(1,bunny_numb),block_number*ones(1,max_sessions_numb),ones(1,2));
CR_crit=0.3;

for i=1:bunny_numb
    sessions_numb=find(cellfun(@(x) isempty(x), CR_amps_ipsi(i,:))==0,1,'last');
    
    for j=1:sessions_numb
        CR_amps_tmp=CR_amps_ipsi{i,j};
        block_size=round(length(CR_amps_tmp)/block_number);
        
        for bl_count=1:block_number
           try 
               CR_amp_block{i,j,1}(bl_count)=nanmean(CR_amps_tmp( 1+(bl_count-1)*block_size:bl_count*block_size,1 ));
               CR_amp_block{i,j,2}(bl_count)=nanmean(CR_amps_tmp( 1+(bl_count-1)*block_size:bl_count*block_size,2 ));
           catch
               CR_amp_block{i,j,1}(bl_count)=nanmean(CR_amps_tmp( 1+(bl_count-1)*block_size:end,1 ));
               CR_amp_block{i,j,2}(bl_count)=nanmean(CR_amps_tmp( 1+(bl_count-1)*block_size:end,2 ));
           end
        end
    end
end
    
figure('units','normalized','outerposition',[0.2 0.4 0.5 0.5]) 
hold on
set(gca,'LineWidth',1)
box off

CR1_amp_conf=reshape(nanstd( cell2mat(CR_amp_block(:,:,1) ) )./sqrt(sum(  isnan(cell2mat(CR_amp_block))==0  )),[],2)';   
CR2_amp_conf=reshape(nanstd( cell2mat(CR_amp_block(:,:,2) ) )./sqrt(sum(  isnan(cell2mat(CR_amp_block))==0  )),[],2)';  

scatter(1/block_number:1/block_number:max_sessions_numb, nanmean(cell2mat(CR_amp_block(:,:,1))), 150, 'k' ,'filled'  )
scatter(1/block_number:1/block_number:max_sessions_numb, nanmean(cell2mat(CR_amp_block(:,:,2))), 150, 'b' ,'filled'  )
errorbar(1/block_number:1/block_number:max_sessions_numb, nanmean(cell2mat(CR_amp_block(:,:,1))),CR1_amp_conf(1,:), CR1_amp_conf(2,:), 'k','MarkerSize',1 ,'LineWidth',1)
errorbar(1/block_number:1/block_number:max_sessions_numb, nanmean(cell2mat(CR_amp_block(:,:,2))),CR2_amp_conf(1,:), CR2_amp_conf(2,:), 'b','MarkerSize',1 ,'LineWidth',1)

xlabel('Session number')
ylabel('CR amplitude, mm')
set(gca,'YTick',0:2:6)
set(gca,'XTick',0:1:3)
axis([0 max_sessions_numb 0 5])
legend(['1st CR'; '2nd CR'])
legend boxoff


% Same analysis for contralateral protocol

figure('units','normalized','outerposition',[0.2 0.4 0.5 0.5]) 
hold on

CR_amps_left_tmp=CR_amps_contra(:,:,1);
CR_amps_left_tmp=cell2mat(CR_amps_left_tmp(:)');
CR_amps_right_tmp=CR_amps_contra(:,:,2);
CR_amps_right_tmp=cell2mat(CR_amps_right_tmp(:)');

scatter(CR_amps_left_tmp,CR_amps_right_tmp,15,'k','filled')

plot([-1 10], [-1 10],'k')
plot([0 10], [0.3 0.3],'k')
plot([0.3 0.3],[0 10],'k')
set(gca,'XTick',0:2:100)
set(gca,'YTick',0:2:100)

axis equal
axis([0 7 0 7])
xlabel('Left CR amplitude, mm')
ylabel('Right CR amplitude, mm')

% conditional probabiltity of right CR amp given left CR amp
amp_bins_L=[0 0.3 2.5:6.5];
amp_bins_R=[0 0.3 2.5:6.5];

prob_CR_amp_right=[];
amp_coord=[];
prob_CR_amp_right_conf=[];

for i=1:length(amp_bins_L)-1
    tr=find(CR_amps_left_tmp>amp_bins_L(i)&CR_amps_left_tmp<amp_bins_L(i+1));
    CR_amp_right_tr=CR_amps_right_tmp(tr);
    prob_CR_amp_right(i,:)=histc(CR_amp_right_tr,amp_bins_R)./sum(histc(CR_amp_right_tr,amp_bins_R));
    amp_coord(i)=(amp_bins_L(i+1)+amp_bins_L(i))/2;
    prob_CR_amp_right_conf(i,:)=bootci(1000,{@(x) (sum(x>0.3))/length(x),CR_amp_right_tr},'alpha',0.05)';
end

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4]) 
hold on

scatter(amp_coord, 100*sum(prob_CR_amp_right(:,2:end),2),100, 'k' ,'filled')
plot(amp_coord, 100*sum(prob_CR_amp_right(:,2:end),2), 'k' )
errorbar( amp_coord, 100*sum(prob_CR_amp_right(:,2:end),2), 100*sum(prob_CR_amp_right(:,2:end),2)-100*prob_CR_amp_right_conf(:,1),100*prob_CR_amp_right_conf(:,2)- 100*sum(prob_CR_amp_right(:,2:end),2), 'k','MarkerSize',1 ,'LineWidth',1)

set(gca,'XTick',0:2:6)
axis([0 6 -1 100])
xlabel('Left CR amplitude, mm')
ylabel('Right CR probability')
set(gca,'YTick',0:20:100)

% Left and Right CR amplitude during Ext. Left Reinf. Right CR sessions 

block_number=10;
bunny_numb=5;
sessions_numb=1;
max_sessions_numb=1;

CR_amp_left=NaN(bunny_numb,max_sessions_numb);
CR_amp_left=mat2cell(CR_amp_left,ones(1,bunny_numb),ones(1,max_sessions_numb));
CR_amp_right=NaN(bunny_numb,max_sessions_numb);
CR_amp_right=mat2cell(CR_amp_right,ones(1,bunny_numb),ones(1,max_sessions_numb));

CR_amp_block_left=[];
CR_amp_block_right=[];

CR_crit_left=0.3;
CR_crit_right=0.3;

for i=1:bunny_numb
    for j=1:sessions_numb
        CR_amp_left=CR_amps_contra{i,j,1};        
        CR_amp_right=CR_amps_contra{i,j,2};              
        block_size=round(length(CR_amp_left)/block_number);
        
        for bl_count=1:block_number
           try 
               CR_amp_block_left{i,j}(bl_count)=mean(CR_amp_left( 1+(bl_count-1)*block_size:bl_count*block_size ));
               CR_amp_block_right{i,j}(bl_count)=mean(CR_amp_right( 1+(bl_count-1)*block_size:bl_count*block_size ));
               CR_perc_block_left{i,j}(bl_count)=100*sum(CR_amp_left( 1+(bl_count-1)*block_size:bl_count*block_size )>CR_crit_left)/block_size;
               CR_perc_block_right{i,j}(bl_count)=100*sum(CR_amp_right( 1+(bl_count-1)*block_size:bl_count*block_size )>CR_crit_right)/block_size;
           catch
               CR_amp_block_left{i,j}(bl_count)=mean(CR_amp_left( 1+(bl_count-1)*block_size:end ));
               CR_amp_block_right{i,j}(bl_count)=mean(CR_amp_right( 1+(bl_count-1)*block_size:end ));
               CR_perc_block_left{i,j}(bl_count)=100*sum(CR_amp_left( 1+(bl_count-1)*block_size:end )>CR_crit_left)/block_size;
               CR_perc_block_right{i,j}(bl_count)=100*sum(CR_amp_right( 1+(bl_count-1)*block_size:end)>CR_crit_right)/block_size;
           end
        end
    end
end


figure('units','normalized','outerposition',[0.2 0.4 0.5 0.5]) 
hold on
set(gca,'LineWidth',1)
box off

CR_amp_left_conf=bootci(1000,{@mean, cell2mat(CR_amp_block_left) },'alpha',0.05);
CR_amp_right_conf=bootci(1000,{@mean, cell2mat(CR_amp_block_right) },'alpha',0.05);

errorbar( 1+block_size/2:block_size:block_number*block_size+block_size/2, mean(cell2mat(CR_amp_block_left)), mean(cell2mat(CR_amp_block_left))-CR_amp_left_conf(1,:), CR_amp_left_conf(2,:) - mean(cell2mat(CR_amp_block_left)), 'k','MarkerSize',1 ,'LineWidth',1)
errorbar( 1+block_size/2:block_size:block_number*block_size+block_size/2, mean(cell2mat(CR_amp_block_right)), mean(cell2mat(CR_amp_block_right))-CR_amp_right_conf(1,:), CR_amp_right_conf(2,:) - mean(cell2mat(CR_amp_block_right)), 'r','MarkerSize',1 ,'LineWidth',1)
scatter(1+block_size/2:block_size:block_number*block_size+block_size/2, mean(cell2mat(CR_amp_block_left)), 150, 'k' ,'filled'  )
scatter(1+block_size/2:block_size:block_number*block_size+block_size/2, mean(cell2mat(CR_amp_block_right)), 150, 'r' ,'filled'  )

xlabel('Trial umber')
ylabel('CR amplitude, mm')
set(gca,'YTick',0:2:6)
set(gca,'XTick',0:25:100)
axis([1 block_number*block_size+block_size/2 -0.1 5])
legend(['Left eye '; 'Right eye'])
legend boxoff

clear

%% Figure 5

% ipsilateral sequence data

load('Figure 5.mat')

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4])
hold on

CS1_CR2_prob_on_CR1s=[];
CS2_CR2_prob_on_CR1s=[];
CS2_CR2_prob_on_nonCR1s=[];

for i=1:length(CS1_CR1_amp_all)
    CS1_CR2_prob_on_CR1s(i)=100*sum(CS1_CR2_amp_all{i}( CS1_CR1_amp_all{i}>=0.3 ) >0.3  )/sum( CS1_CR1_amp_all{i}>=0.3 );
    CS2_CR2_prob_on_CR1s(i)=100*sum(CS2_CR2_amp_all{i}( CS2_CR1_amp_all{i}>=0.3) >0.3  )/sum( CS2_CR1_amp_all{i}>=0.3 );
    CS2_CR2_prob_on_nonCR1s(i)=100*sum(CS2_CR2_amp_all{i}( CS2_CR1_amp_all{i}<0.3 ) >0.3  )/sum( CS2_CR1_amp_all{i}<0.3 );
end

[r, p]=corrcoef(CS1_CR2_prob_on_CR1s,CS2_CR2_prob_on_CR1s,'rows','complete');

scatter(CS1_CR2_prob_on_CR1s,CS2_CR2_prob_on_CR1s,50,'filled','k')

plot([0 100],[0 100],'color',[0.5 0.5 0.5])
axis square
axis([0 100 0 100])
set(gca,'XTick',0:25:100)
set(gca,'YTick',0:25:100)
xlabel('CS1, 2nd CR probability')
ylabel('CS2, 2nd CR probability')


CS1_CR2_prob_on_CR1s_conf=bootci(1000,{@nanmean, CS1_CR2_prob_on_CR1s},'alpha',0.05)';
CS2_CR2_prob_on_CR1s_conf=bootci(1000,{@nanmean, CS2_CR2_prob_on_CR1s},'alpha',0.05)';
CS2_CR2_prob_on_nonCR1s_conf=bootci(1000,{@nanmean, CS2_CR2_prob_on_nonCR1s},'alpha',0.05)';

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4])
hold on
set(gca,'LineWidth',1)
box off

bar(1, nanmean(CS1_CR2_prob_on_CR1s),0.5,'FaceColor',[1 1 1])
scatter(ones(1,length(CS1_CR2_prob_on_CR1s)),CS1_CR2_prob_on_CR1s,20, 'b' )
errorbar( 1, nanmean(CS1_CR2_prob_on_CR1s), nanmean(CS1_CR2_prob_on_CR1s)-CS1_CR2_prob_on_CR1s_conf(1),CS1_CR2_prob_on_CR1s_conf(2)- nanmean(CS1_CR2_prob_on_CR1s), 'b','MarkerSize',1 ,'LineWidth',1)

bar(2, nanmean(CS2_CR2_prob_on_CR1s),0.5,'FaceColor',[1 1 1])
scatter(2*ones(1,length(CS2_CR2_prob_on_CR1s)),CS2_CR2_prob_on_CR1s,20, [0.2 0 0])
errorbar( 2, nanmean(CS2_CR2_prob_on_CR1s), nanmean(CS2_CR2_prob_on_CR1s)-CS2_CR2_prob_on_CR1s_conf(1),CS2_CR2_prob_on_CR1s_conf(2)- nanmean(CS2_CR2_prob_on_CR1s), 'color', [0.2 0 0],'MarkerSize',1 ,'LineWidth',1)

bar(3, nanmean(CS2_CR2_prob_on_nonCR1s),0.5,'FaceColor',[1 1 1])
scatter(3*ones(1,length(CS2_CR2_prob_on_nonCR1s)),CS2_CR2_prob_on_nonCR1s,20, 'k' )
errorbar( 3, nanmean(CS2_CR2_prob_on_nonCR1s), nanmean(CS2_CR2_prob_on_nonCR1s)-CS2_CR2_prob_on_nonCR1s_conf(1),CS2_CR2_prob_on_nonCR1s_conf(2)- nanmean(CS2_CR2_prob_on_nonCR1s), 'k','MarkerSize',1 ,'LineWidth',1)

axis([0 4 0 100])
ylabel('2nd CR probability')
set(gca,'YTick',0:20:100)
set(gca,'XTick',0:1:10)

% contralateral sequence data

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4])
hold on

CS1_R_CR_prob_on_L_CRs=[];
CS2_R_CR_prob_on_L_CRs=[];
CS2_R_CR_prob_on_L_nonCRs=[];

for i=1:length(CS1_CR_amp_left_all)
    CS1_R_CR_prob_on_L_CRs(i)=100*sum(CS1_CR_amp_right_all{i}( CS1_CR_amp_left_all{i}>3 ) >0.3  )/sum( CS1_CR_amp_left_all{i}>3 );
    CS2_R_CR_prob_on_L_CRs(i)=100*sum(CS2_CR_amp_right_all{i}( CS2_CR_amp_left_all{i}>3 ) >0.3  )/sum( CS2_CR_amp_left_all{i}>3 );
    CS2_R_CR_prob_on_L_nonCRs(i)=100*sum(CS2_CR_amp_right_all{i}( CS2_CR_amp_left_all{i}<0.3 ) >0.3  )/sum( CS2_CR_amp_left_all{i}<0.3 );
end

[r, p]=corrcoef(CS1_R_CR_prob_on_L_CRs,CS2_R_CR_prob_on_L_CRs,'rows','complete');
scatter(CS1_R_CR_prob_on_L_CRs,CS2_R_CR_prob_on_L_CRs,50,'filled','k')
plot([0 100],[0 100],'color',[0.5 0.5 0.5])
axis square
axis([0 100 0 100])
set(gca,'XTick',0:25:100)
set(gca,'YTick',0:25:100)
xlabel('CS1, Right CR probability')
ylabel('CS2, Right CR probability')


figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4])
hold on
set(gca,'LineWidth',1)
box off

CS1_R_CR_prob_on_L_CRs_conf=bootci(2000,{@nanmean, CS1_R_CR_prob_on_L_CRs},'alpha',0.05)';
CS2_R_CR_prob_on_L_CRs_conf=bootci(2000,{@nanmean, CS2_R_CR_prob_on_L_CRs},'alpha',0.05)';
CS2_R_CR_prob_on_L_nonCRs_conf=bootci(2000,{@nanmean, CS2_R_CR_prob_on_L_nonCRs},'alpha',0.05)';

bar(1, nanmean(CS1_R_CR_prob_on_L_CRs),0.5,'FaceColor',[1 1 1])
scatter(ones(1,length(CS1_R_CR_prob_on_L_CRs)),CS1_R_CR_prob_on_L_CRs,20, 'r' )
errorbar( 1, nanmean(CS1_R_CR_prob_on_L_CRs), nanmean(CS1_R_CR_prob_on_L_CRs)-CS1_R_CR_prob_on_L_CRs_conf(1),CS1_R_CR_prob_on_L_CRs_conf(2)- nanmean(CS1_R_CR_prob_on_L_CRs), 'r','MarkerSize',1 ,'LineWidth',1)

bar(2, nanmean(CS2_R_CR_prob_on_L_CRs),0.5,'FaceColor',[1 1 1])
scatter(2*ones(1,length(CS2_R_CR_prob_on_L_CRs)),CS2_R_CR_prob_on_L_CRs,20, [0.2 0 0])
errorbar( 2, nanmean(CS2_R_CR_prob_on_L_CRs), nanmean(CS2_R_CR_prob_on_L_CRs)-CS2_R_CR_prob_on_L_CRs_conf(1),CS2_R_CR_prob_on_L_CRs_conf(2)- nanmean(CS2_R_CR_prob_on_L_CRs), 'color', [0.2 0 0],'MarkerSize',1 ,'LineWidth',1)

bar(3, nanmean(CS2_R_CR_prob_on_L_nonCRs),0.5,'FaceColor',[1 1 1])
scatter(3*ones(1,length(CS2_R_CR_prob_on_L_nonCRs)),CS2_R_CR_prob_on_L_nonCRs,20, 'k' )
errorbar( 3, nanmean(CS2_R_CR_prob_on_L_nonCRs), nanmean(CS2_R_CR_prob_on_L_nonCRs)-CS2_R_CR_prob_on_L_nonCRs_conf(1),CS2_R_CR_prob_on_L_nonCRs_conf(2)- nanmean(CS2_R_CR_prob_on_L_nonCRs), 'k','MarkerSize',1 ,'LineWidth',1)

axis([0 4 0 100])
ylabel('Right CR probability')
set(gca,'YTick',0:20:100)
set(gca,'XTick',0:1:10)

title('CS1_CR1s CS2_CR1s CS2_nonCR1s','interpreter','none')

clear

%% Figure 6

load('Figure 6.mat')

% ipsilateral sequence data

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4]) 
hold on
set(gca,'LineWidth',1)
box off

first_CR_amp=CR_amps_ipsi(:,1);
second_CR_amp=CR_amps_ipsi(:,2);
scatter(CR_amps_ipsi(:,1),CR_amps_ipsi(:,2),25,'k','filled')

plot([-1 10], [-1 10],'k')
plot([0 10], [0.3 0.3],'k')
plot([0.3 0.3],[0 10],'k')
axis equal
axis([-0. 7 -0. 7])
set(gca,'YTick',0:2:6)
set(gca,'XTick',0:2:6)
xlabel('First CR amplitude, mm')
ylabel('Second CR amplitude, mm')

%conditional probabiltity of 2nd CR  given 1st CR amp

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4]) 
hold on
set(gca,'LineWidth',1)
box off

amp_coord=[];
prob_2nd_CR_conf=[];
prob_2nd_CR=[];

amp_bins=[0 0.3 2.5:6.5];

for i=1:length(amp_bins)-1
    p_count=0;
    tr=find(first_CR_amp>amp_bins(i)&first_CR_amp<amp_bins(i+1));
    second_CR_amp_tr=second_CR_amp(tr);
    prob_2nd_CR(i)=sum(second_CR_amp_tr>0.3)/length(second_CR_amp_tr);
    amp_coord(i)=(amp_bins(i+1)+amp_bins(i))/2;
    prob_2nd_CR_conf(i,:)=bootci(1000,{@(x) (sum(x>0.3))/length(x),second_CR_amp_tr},'alpha',0.05);
end

scatter(amp_coord, 100*prob_2nd_CR,100, 'k' ,'filled')
plot(amp_coord, 100*prob_2nd_CR, 'k' )
errorbar( amp_coord, 100*prob_2nd_CR, 100*prob_2nd_CR-100*prob_2nd_CR_conf(:,1)',100*prob_2nd_CR_conf(:,2)'- 100*prob_2nd_CR, 'k','MarkerSize',1 ,'LineWidth',1)

axis([0 6 0 100])
xlabel('First CR amplitude, mm')
ylabel('Second CR probability')
set(gca,'YTick',0:20:100)

% contralateral sequence data

figure('units','normalized','outerposition',[0.2 0.4 0.5 0.5]) 
set(gca,'LineWidth',1)
box off
hold on
scatter(comb_CR_amp_contra(1,:), comb_CR_amp_contra(2,:),25, 'k' ,'filled')
plot([-1 10], [-1 10],'k')
plot([-1 10], [0.3 0.3],'k')
plot([0.3 0.3],[-1 10],'k')
axis equal
axis([-0 7 -0 7])
xlabel('Left CR amplitude, mm')
ylabel('Right CR amplitude, mm')
set(gca,'YTick',0:2:10)
set(gca,'XTick',0:2:10)

amp_bins_L=[0 0.3 2.5:6.5];
amp_bins_R=[0 0.3 6];
prob_CR_amp_right=[];
prob_CR_amp_right_conf=[];
amp_coord=[];

for i=1:length(amp_bins_L)-1
    
    tr=find(comb_CR_amp_contra(1,:)>amp_bins_L(i)&comb_CR_amp_contra(1,:)<amp_bins_L(i+1));
    CR_amp_right_tr=comb_CR_amp_contra(2,tr);
    prob_CR_amp_right(i,:)=histc(CR_amp_right_tr,amp_bins_R)./sum(histc(CR_amp_right_tr,amp_bins_R));
    amp_coord(i)=(amp_bins_L(i+1)+amp_bins_L(i))/2;
    prob_CR_amp_right_conf(i,:)=bootci(2000,{@(x) (sum(x>0.3))/length(x),CR_amp_right_tr},'alpha',0.05)';

end

figure('units','normalized','outerposition',[0.2 0.4 0.4 0.4])
hold on
set(gca,'LineWidth',1)
box off
scatter(amp_coord, 100*sum(prob_CR_amp_right(:,2:end),2),100, 'k' ,'filled')
errorbar( amp_coord, 100*sum(prob_CR_amp_right(:,2:end),2), 100*sum(prob_CR_amp_right(:,2:end),2)-100*prob_CR_amp_right_conf(:,1),100*prob_CR_amp_right_conf(:,2)- 100*sum(prob_CR_amp_right(:,2:end),2), 'k','MarkerSize',1 ,'LineWidth',1)
axis([0 6 0 100])
xlabel('Left CR amplitude, mm')
ylabel('Right CR probability')
set(gca,'YTick',0:20:100)



%% Figure 8

load('Figure 8.mat')

% PC firing rate on nonCRs, early, late trials: 1st CRs

ISI1=500; 
delta=1500; %in ms, baseline before CS

first_CR_onset_time=[first_US_trials.CR_onset second_US_trials.CR_onset CS_alone.CR_onset ];
first_CR_amp=[(first_US_trials.eyelid_pos(200+ISI1,:)) (second_US_trials.eyelid_pos(200+ISI1,:)) (CS_alone.eyelid_pos(200+ISI1,:)) ];

eyelid_pos_tmp=[first_US_trials.eyelid_pos second_US_trials.eyelid_pos CS_alone.eyelid_pos];
fr_tmp=[first_US_trials.eyelid_PC_firing_rate ; second_US_trials.eyelid_PC_firing_rate ; CS_alone.eyelid_PC_firing_rate];

nonCR_trials=find(first_CR_amp<0.3);

[sorted_onset,ind]=sort(first_CR_onset_time);
sorted_onset(isnan(sorted_onset)==1)=[];
mean_1st_CR_onset=sorted_onset(round(length(sorted_onset)/2));
group_1_CRs=find(first_CR_onset_time<mean_1st_CR_onset);
group_2_CRs=find(first_CR_onset_time>=mean_1st_CR_onset);

figure('units','normalized','outerposition',[0.1 0.1 0.25 0.7]) 
subplot(2,1,1)
hold on
    
plot(-199:ISI1,mean(eyelid_pos_tmp(1:200+ISI1,nonCR_trials),2),'k','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI1,nonCR_trials)');
plot(-199:ISI1,muci(1,1:ISI1+200),'k','LineWidth',1)
plot(-199:ISI1,muci(2,1:ISI1+200),'k','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-199:ISI1,mean(eyelid_pos_tmp(1:200+ISI1,group_1_CRs),2),'r','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI1,group_1_CRs)');
plot(-199:ISI1,muci(1,1:ISI1+200),'r','LineWidth',1)
plot(-199:ISI1,muci(2,1:ISI1+200),'r','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-199:ISI1,mean(eyelid_pos_tmp(1:200+ISI1,group_2_CRs),2),'b','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI1,group_2_CRs)');
plot(-199:ISI1,muci(1,1:ISI1+200),'b','LineWidth',1)
plot(-199:ISI1,muci(2,1:ISI1+200),'b','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot([0 0],[-1 10],'k')
plot([ISI1 ISI1],[-1 10],'k')

set(gca,'LineWidth',1)
set(gca,'XTick',-500:500:2000)
set(gca,'YTick',0:2:8)
axis([-300 ISI1 -0.2 6])

subplot(2,1,2)
hold on

plot(-delta+1:ISI1,mean(fr_tmp(nonCR_trials,1:delta+ISI1)),'k','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(nonCR_trials,:));
plot(-delta+1:ISI1,muci(1,1:ISI1+delta),'k','LineWidth',1)
plot(-delta+1:ISI1,muci(2,1:ISI1+delta),'k','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-delta+1:ISI1,mean(fr_tmp(group_1_CRs,1:delta+ISI1)),'r','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(group_1_CRs,:));
plot(-delta+1:ISI1,muci(1,1:ISI1+delta),'r','LineWidth',1)
plot(-delta+1:ISI1,muci(2,1:ISI1+delta),'r','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-delta+1:ISI1,mean(fr_tmp(group_2_CRs,1:delta+ISI1)),'b','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(group_2_CRs,:));
plot(-delta+1:ISI1,muci(1,1:ISI1+delta),'b','LineWidth',1)
plot(-delta+1:ISI1,muci(2,1:ISI1+delta),'b','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot([0 0],[0 2],'k')
plot([ISI1 ISI1],[0 2],'k')
plot([-500 ISI1],[1 1],'k')
set(gca,'LineWidth',1)
set(gca,'XTick',-500:500:2000)
set(gca,'YTick',0:0.2:4)
axis([-300 ISI1 0.3 1.2])


% PC firing on nonCRs, early, late trials: 2nd CRs

ISI2=1100;

eyelid_pos_tmp=[second_US_trials.eyelid_pos CS_alone.eyelid_pos ];
fr_tmp=[second_US_trials.eyelid_PC_firing_rate ; CS_alone.eyelid_PC_firing_rate];

second_CR_onset=[];

for q=1:length(eyelid_pos_tmp(1,:))
    pos_tr=eyelid_pos_tmp(1:1301,q);
    pos_tr(end)=0;
    [pks,locs] = findpeaks(pos_tr,'MinPeakProminence',0.5,'MinPeakDistance',200);
    if length(pks)==2&&locs(2)>(ISI2+200-100)
        second_CR_amp(q)=pks(2);
        [~,locs] = findpeaks(-pos_tr,'MinPeakProminence',0.5,'MinPeakDistance',200);
        try
            second_CR_onset(q)=locs(1)-200;
        catch
            second_CR_onset(q)=NaN;
        end
    else
        if length(pks)==1
            second_CR_amp(q)=0;
        else
            second_CR_amp(q)=NaN;
        end
        second_CR_onset(q)=NaN;
    end
end

[sorted_onset,ind]=sort(second_CR_onset);
sorted_onset(isnan(sorted_onset)==1)=[];
mean_2nd_CR_onset=sorted_onset(round(length(sorted_onset)/2));
group_1_CRs=find(second_CR_onset<mean_2nd_CR_onset);
group_2_CRs=find(second_CR_onset>=mean_2nd_CR_onset);

nonCR_trials=find(second_CR_amp<0.3);

figure('units','normalized','outerposition',[0.1 0.1 0.5 0.7]) 
subplot(2,1,1)
hold on
    
plot(-199:ISI2,mean(eyelid_pos_tmp(1:200+ISI2,nonCR_trials),2),'k','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI2,nonCR_trials)');
plot(-199:ISI2,muci(1,1:ISI2+200),'k','LineWidth',1)
plot(-199:ISI2,muci(2,1:ISI2+200),'k','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-199:ISI2,mean(eyelid_pos_tmp(1:200+ISI2,group_1_CRs),2),'r','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI2,group_1_CRs)');
plot(-199:ISI2,muci(1,1:ISI2+200),'r','LineWidth',1)
plot(-199:ISI2,muci(2,1:ISI2+200),'r','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-199:ISI2,mean(eyelid_pos_tmp(1:200+ISI2,group_2_CRs),2),'b','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(eyelid_pos_tmp(1:200+ISI2,group_2_CRs)');
plot(-199:ISI2,muci(1,1:ISI2+200),'b','LineWidth',1)
plot(-199:ISI2,muci(2,1:ISI2+200),'b','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot([0 0],[-1 10],'k')
plot([ISI1 ISI1],[-1 10],'k')

set(gca,'LineWidth',1)
set(gca,'XTick',-500:500:2000)
set(gca,'YTick',0:2:8)
axis([-300 ISI2 -0.2 6])

subplot(2,1,2)
hold on

plot(-delta+1:ISI2,mean(fr_tmp(nonCR_trials,1:delta+ISI2)),'k','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(nonCR_trials,:));
plot(-delta+1:ISI2,muci(1,1:ISI2+delta),'k','LineWidth',1)
plot(-delta+1:ISI2,muci(2,1:ISI2+delta),'k','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-delta+1:ISI2,mean(fr_tmp(group_1_CRs,1:delta+ISI2)),'r','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(group_1_CRs,:));
plot(-delta+1:ISI2,muci(1,1:ISI2+delta),'r','LineWidth',1)
plot(-delta+1:ISI2,muci(2,1:ISI2+delta),'r','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot(-delta+1:ISI2,mean(fr_tmp(group_2_CRs,1:delta+ISI2)),'b','LineWidth',2)
[muhat,sigmahat,muci,sigmaci] = normfit(fr_tmp(group_2_CRs,:));
plot(-delta+1:ISI2,muci(1,1:ISI2+delta),'b','LineWidth',1)
plot(-delta+1:ISI2,muci(2,1:ISI2+delta),'b','LineWidth',1)
clearvars muhat sigmahat muci sigmaci

plot([0 0],[0 2],'k')
plot([ISI1 ISI1],[0 2],'k')
plot([-500 ISI2],[1 1],'k')

set(gca,'LineWidth',1)
set(gca,'XTick',-500:500:2000)
set(gca,'YTick',0:0.2:4)
axis([-300 ISI2 0.3 1.2])

clear

























