clc
clear all
close all
%%
%inputs
%CHANGES TO MAKE-values of field length are in voxels ,have to convert
vl=0.25;%voxel length
%firr_thresh=0.2;%0.2
tilted=1;%check if tilted
%calculate stats from which neuron to which neuron?
if tilted
%load("Tilted_data_props_b0.9b_std2A.mat")
load("Tilted_data_props_b0.8pi_std3d1.5_std2d1.5.mat")
%load("field_elongstats_tilted.mat")
%neuron_end = 50;
neuron_start = 1;
neurons = [5,6,7,10,14,18,25,26,33,36,39,41,46];
%neurons = [2 6 7 10 11 12 19 21 24 25 27 30 32 33 37 41 45 47 48 50]; % with 0.9 std threshold
%neurons = [0, 1, 2, 3, 5, 6, 7, 8, 10, 11, 14, 19, 20, 23, 24, 25, 26, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 40, 41, 42, 43, 45, 46, 47, 48, 49]+1; %With aligned phi case
%neurons = [ 1 3 4 6 9 10 11 13 14 15 17 21 26 27 32 33 34 35 37 39 41 42 43 44 46 48 49];
else
load("Aligned_data_props_b2.51_t1.7_2Dt1.5.mat")
neuron_start = 1;
%load("field_elongstats_aligned.mat")
%neuron_end = 47;
neurons = [2, 4, 6, 8, 11, 16, 18, 24, 26, 30, 33, 34, 41] + 1;
end
minv=50;%minimum voxels to be considered as place field
%CHOOSE PLOTS TO PLOT
identify_placef=1;%set this to 1 to identify place fields ,check min voxels
calc_orientation_elongation=1;
plot_orientation=1; %plot orientation on sphere
plot_elongation=1;
plot_centroids=1;
plot_binary=1;%binary morphology
fields_per_cell=0;%plot fields per cell
plot_pf=0;%plot place field firing rate map
calc_proportion=0; %calculate fields parallel to x,y,z
plot_flength=0;
plotfirr=0;
neuron_firr=1;%neuron for which to plot firr
rotm = [[ 0.57357644, 0., 0.81915204],
[ 0.57922797, 0.70710678, -0.40557979],
[-0.57922797, 0.70710678, 0.40557979]];
% rotm = [[ 0.70710678, -0.57922797, -0.40557979],
% [ 0. , 0.57357644, -0.81915204],
% [ 0.70710678, 0.57922797, 0.40557979]];
XYZ = [1 0 0; 0 1 0; 0 0 1];
ABC = rotm*XYZ;
axes_cords = [XYZ; ABC];
% set(gca,'FontSize',30);
%% %%%%%%% IDENTIFY PLACE CELLS %%%%%%%%%%%%
x1=1.125:0.25:5.875;
y1=1.125:0.25:5.875;
z1=1.125:0.25:5.875;
[X,Y,Z]=meshgrid(x1, y1 ,z1);
voxel_coords=[Y(:) X(:) Z(:)];%voxel coords increase along x then y then z
orientation=double.empty;
orientationplot=double.empty;
elongation=double.empty;
yfl=double.empty;
xfl=double.empty;
zfl=double.empty;
centroids=double.empty;
%%identify place fields from all neurons and store in pf structure
if identify_placef
pf=struct;
pf_counter=1;
% for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
for neuron = neurons
propstable=FIRR_3dprops(neuron).props; %get 3d properties table of the neuron
firrtable=FIRR_3dprops(neuron).firr;
for i=1:height(propstable)
vol=propstable{i,"Volume"};
if vol>minv %if volume of cluster is greater than minvoxels
pf(pf_counter).props=propstable(i,:);
pf(pf_counter).neuron=neuron;
pf(pf_counter).firr=firrtable;
pf_counter=pf_counter+1;
end
end
end
end
%% %%%%%%% PLOT CENTROIDS FIG 4 %%%%%%%%%%%%
if plot_centroids
%%%%%%%%%%%%% FIG 4 a1) Aligned or FIG 4 b1) Tilted %%%%%%%%%%%%%%
set(0, 'DefaultFigureRenderer', 'painters');
for i=1:size(pf,2)
centroid=pf(i).props.Centroid;
centroids=[centroids;centroid];
centroids=centroids;
end
% centroids = rotm * centroids';
centroids_scaled=(5/20)*(centroids)+1;
if tilted == 1
centroids_scaled1 = rotm * centroids_scaled';
centroids_scaled1 = centroids_scaled1';
% x_centroid = [10 0 0];
% x_centroid_rot = x_centroid * rotm;
% y_centroid = [0 10 0];
% y_centroid_rot = y_centroid * rotm;
% z_centroid = [0 0 10];
% z_centroid_rot = z_centroid * rotm;
end
%%hold off
%figure
if tilted == 1
V1 = [1 1 1];
V2 = [6 1 1];
V3 = [6 1 6];
V4 = [1 1 6];
V5 = [1 6 6];
V6 = [1 6 1];
V7 = [1 6 6];
V8 = [6 6 6];
V9 = [6 6 1];
V10 = [1 6 1];
V11 = [1 1 1];
V12 = [1 1 6];
V13 = [6 1 6];
V14 = [6 6 6];
V15 = [6 6 1];
V16 = [6 1 1];
C = [V1;V2;V3;V4;V5;V6;V7;V8;V9;V10;V11;V12;V13;V14;V15;V16];
E = rotm*C';
E = E';
figure
plot3(E(:,1),E(:,2),E(:,3),'black')
hold on
end
if tilted
scatter3(centroids_scaled1(:,1) ,centroids_scaled1(:,2), centroids_scaled1(:,3),30, centroids_scaled1(:,3),"filled")
title("Centroids of place fields- Tilted lattice")
centroid = mean(E);
% Define line directions in ABC axes
line_directions = [ABC(:, 1)'; ABC(:, 2)'; ABC(:, 3)'];
% Define line length
line_length = 5;
% Calculate line end points
line_ends = centroid + line_length * line_directions;
% Plot the lines
hold on;
line_colors = {'red', 'green', 'blue'};
for i = 1:3
line_coords = [centroid; line_ends(i, :)];
plot3(line_coords(:, 1), line_coords(:, 2), line_coords(:, 3), 'Color', line_colors{i},'LineWidth', 2 );
end
% hold on
% plot3(x_centroid_rot, y_centroid_rot, z_centroid_rot);
daspect([1 1 1]);
%hold on
else
scatter3(centroids_scaled(:,1) ,centroids_scaled(:,2), centroids_scaled(:,3),30, centroids_scaled(:,3),"filled")
title("Centroids of place fields- Aligned lattice")
%hold on
end
xlabel('x','Fontsize',15)
ylabel('y','Fontsize',15)
zlabel('z','Fontsize',15)
if tilted
medx=mean(centroids_scaled1(:,1));
medy=mean(centroids_scaled1(:,2));
medz=mean(centroids_scaled1(:,3));
else
medx=mean(centroids_scaled(:,1));
medy=mean(centroids_scaled(:,2));
medz=mean(centroids_scaled(:,3));
end
if tilted == 0
hold on
%
scatter3(medx,medy,medz,100,'filled','red')
hold on
end
if tilted == 1
hold on
scatter3(medx,medy,medz,100,'filled','red')% scatter3(medx,3.5,3.5,100,'filled','green')
end
x2=3.5*ones(length(x1),1);% x2=3.5*ones(length(x1),1);
y2=3.5*ones(length(y1),1);
z2=3.5*ones(length(z1),1);
xyz_rot = rotm * [x2 y2 z2]';
if tilted ==1
% plot3(x1, xyz_rot(2, :), xyz_rot(3, :), 'Color', 'red');
% hold on
% plot3(xyz_rot(1, :), y1, xyz_rot(3, :), 'Color', 'blue');
% hold on
% plot3(xyz_rot(1, :), xyz_rot(2, :), z1, 'Color', 'green');
% hold on
% plot3(x1,y2,z2,'Color','black');
% hold on
% plot3(x2,y1,z2,'Color','black');
% hold on
% plot3(x2,y2,z1,'Color','black');
% hold on
end
hold off
%%%%%%%%%%% Fig 4a2) Aligned or Fig 4b2) Tilted %%%%%%%%%%%
if tilted
cc = inv(rotm)*centroids_scaled';
cc = cc';
end
cc = centroids_scaled;
center = [3.5, 3.5, 3.5];
cc = cc - center;
bw = 0.7;
figure;
xi = -5:0.1:5;
[~,~,~,f1] = histDENSITY(cc(:,1),xi,bw);
[~,~,~,f2] = histDENSITY(cc(:,2),xi,bw);
[~,~,~,f3] = histDENSITY(cc(:,3),xi,bw);
leng_line = max([max(f1) max(f2) max(f3)]);
p1 = plot(xi(:),f1,'LineWidth',1.5, 'color', 'red'); hold on;
p2 = plot(xi(:),f2,'LineWidth',1.5, 'Color', 'green');
p3 = plot(xi(:),f3,'LineWidth',1.5, 'Color', 'blue');
line([nanmedian(cc(:,1)),nanmedian(cc(:,1))], [0 leng_line], 'color', 'red');
line([nanmedian(cc(:,2)),nanmedian(cc(:,2))], [0 leng_line], 'color', 'green');
line([nanmedian(cc(:,3)),nanmedian(cc(:,3))], [0 leng_line], 'color', 'blue');
legend('A', 'B', 'C');
xlabel('Distance from center'); ylabel('Proportion of fields');
%%%%%%% Fig 4a3) Aligned or Fig 4b3) Tilted %%%%%%%%
xi = -0.5:.001:0.5;
med_rand_x =[]; med_rand_y =[]; med_rand_z=[];
for i = 1:1000
med_rand_x = [med_rand_x nanmedian((5*rand(1,length(cc)))-2.5)];
med_rand_y = [med_rand_y nanmedian((5*rand(1,length(cc)))-2.5)];
med_rand_z = [med_rand_z nanmedian((5*rand(1,length(cc)))-2.5)];
end
med_rand = [med_rand_x; med_rand_y; med_rand_z]';
colors = ['r'; 'g'; 'b'];
figure;
for j =1:3
subplot(3,1,j);
temp = med_rand(:,j);
[~,~,~,f] = histDENSITY(temp,xi,0.01);
p1 = plot(xi(:),f,'LineWidth',1.5, 'color', colors(j)); %hold on;
line([nanmedian(cc(:,j)),nanmedian(cc(:,j))], [0 max(f)], 'color', 'black', 'LineWidth', 1.5);
line([prctile(temp, 97.5), prctile(temp, 97.5)], [0 max(f)], 'color', 'red', 'LineWidth', 1.5);
line([prctile(temp, 2.5), prctile(temp, 2.5)], [0 max(f)], 'color', 'red','LineWidth', 1.5);
xlabel("Median field position");
ylabel("Proportion of fields");
end
end
%% %%%%%%% CALCULATE ORIENTATION AND ELONGATION %%%%%%%%%%%%
%orientation ,elongation
if calc_orientation_elongation
orientation=double.empty;
for i=1:length(pf)%calculate stats
% i
%calculate orientation
b=pf(i).props.ConvexHull{1,1};
a=pca(b);
if tilted
a = rotm * a;
orientation=[orientation a(:,1)];
orientationplot=[orientationplot a(:,1) -a(:,1)];
else
orientation=[orientation a(:,1)];
orientationplot=[orientationplot a(:,1) -a(:,1)];
end
%calculate elongation
t=pf(i).props;
arr=t.PrincipalAxisLength;
elong=2*arr(1)/(arr(2)+arr(3));
elongation=[elongation ;elong];
end
end
extent_field = [];
%for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
for neuron = neurons
propstable=FIRR_3dprops(neuron).props;
for i=1:height(propstable)
vol=propstable{i,"Volume"};
if vol>minv
extent_field = [extent_field; propstable.BoundingBox(i,4:6)];
end
end
end
figure
if tilted
xi = 0:1:25;
window_size = 10;
[~,~,~,f1] = histDENSITY(extent_field(:,1),xi,1);
[~,~,~,f2] = histDENSITY(extent_field(:,2),xi,1);
[~,~,~,f3] = histDENSITY(extent_field(:,3),xi,1);
smoothed_f1 = movmean(f1, window_size);
smoothed_f2 = movmean(f2, window_size);
smoothed_f3 = movmean(f3, window_size);
p1 = plot(xi(:), smoothed_f1, 'color', 'r', 'LineWidth', 1.5);hold on;
p2 = plot(xi(:), smoothed_f2, 'color', 'g', 'LineWidth', 1.5);
p3 = plot(xi(:), smoothed_f3, 'color', 'b', 'LineWidth', 1.5);
xlabel("Field extent");
ylabel("Fields (%)");
legend('A', 'B', 'C');
else
xi = 0:1:25;
[~,~,~,f1] = histDENSITY(extent_field(:,1),xi,2.5);
[~,~,~,f2] = histDENSITY(extent_field(:,2),xi,2.5);
[~,~,~,f3] = histDENSITY(extent_field(:,3),xi,2.5);
p1 = plot(xi(:), f1, 'color', 'r', 'LineWidth', 1.5);hold on;
p2 = plot(xi(:), f2, 'color', 'g', 'LineWidth', 1.5);
p3 = plot(xi(:), f3, 'color', 'b', 'LineWidth', 1.5);
end
%% %%%%%%% BINARY MORPHOLOGY FIG 7 %%%%%%%%%%%%
if plot_binary
total=0;
seq_len = 19;
%create table like array that has SE length and correponding
%proportions
Px=double.empty;
Py=double.empty;
Pz=double.empty;
% Px(1,:)=[3:2:19];
% Py(1,:)=[3:2:19];
% Pz(1,:)=[3:2:19];
Px(1,:)=1:1:seq_len;
Py(1,:)=1:1:seq_len;
Pz(1,:)=1:1:seq_len;
Px(2,:)=zeros(1,seq_len);
Py(2,:)=zeros(1,seq_len);
Pz(2,:)=zeros(1,seq_len);
%seq_len_ratio = [];
Seq_len_ratio = [];
ind=1 ; %index simply represents index of array l which contains structuring element length ;
for len=1:1:seq_len %
ex = strel('cuboid',[len 1 1]); %elementx
ey = strel('cuboid',[1 len 1]);
ez = strel('cuboid',[1 1 len]);
total=0;
%for neuron=setdiff(neuron_start:neuron_end,avoid_neurons) %binarize rate map of each neuron using 10% of maximum as thresh and make binary volumetric image
%neuron
for neuron = neurons
bw=zeros(20,20,20);
firr=FIRR_3dprops(neuron).firr;
counter=1;
threshold=0.1*max(firr);
for k=1:20 % this loop is to create binary volumetric image for a neuron (with threshold)
j=1;
for j=1:20
i=1;
for i=1:20
%counter
if firr(counter)>threshold
bw(i,j,k)=1;
end
counter=counter+1;
end
end
end
Bx=imerode(bw,ex);
By=imerode(bw,ey);
Bz=imerode(bw,ez);
total1=sum(Bx,"all")+sum(By,"all")+sum(Bz,"all");%total number of voxels with a particular connectivity for a particular neuron's firing rate
total=total+total1; %total ,after complete iteration of all neurons,is the sum of all voxels remaining for that particular SE length
if neuron==neuron_start
Px(2,ind)=sum(Bx,"all");
Py(2,ind)=sum(By,"all");
Pz(2,ind)=sum(Bz,"all");
else
Px(2,ind)=Px(2,ind)+sum(Bx,"all");%sum of remaining voxels of every iterated neuron along x axis
Py(2,ind)=Py(2,ind)+sum(By,"all");
Pz(2,ind)=Pz(2,ind)+sum(Bz,"all");
end
neurons_ratio{neuron,1} = len;
neurons_ratio{neuron,2} = sum(Bx,"all")/total1;
neurons_ratio{neuron,3} = sum(By,"all")/total1;
neurons_ratio{neuron,4} = sum(Bz,"all")/total1;
end%neuron is inner loop
% if total == 0
% total = 0.00000001;
% end
% total = total + 0.0000001;
Px(2,ind)=Px(2,ind)/total;
Px
Py(2,ind)=Py(2,ind)/total;
Pz(2,ind)=Pz(2,ind)/total;
ind=ind+1;
%neurons_ratio = neurons_ratio(~isnan(neurons_ratio));
Seq_len_ratio = vertcat(Seq_len_ratio, neurons_ratio);
%seq_len_ratio{len} = neurons_ratio;
end %strel is outer loop
Seq_len_ratio = cell2mat(Seq_len_ratio);
%Seq_len_ratio = Seq_len_ratio(~isnan(Seq_len_ratio(:,2)), :);
%Seq_len_ratio = Seq_len_ratio(~cellfun(@isempty, Seq_len_ratio), :);
%Seq_len_ratio = Seq_len_ratio(~any(isnan(Seq_len_ratio), 2), :);
figure;
plot(Px(1,:),Px(2,:), 'Color', 'r', 'LineWidth', 1.5)
hold on
plot(Py(1,:),Py(2,:), 'Color', 'g', 'LineWidth', 1.5)
hold on
plot(Pz(1,:),Pz(2,:), 'Color', 'b', 'LineWidth', 1.5)
ylabel('Proportion of voxels','FontSize',15)
xlabel('Structuring element length','FontSize',15)
if tilted
legend('A','B','C')
else
legend('X','Y','Z')
end
xlim([0 20]); ylim([0 1]);
title('Binary morphology','Fontsize',16)
end
%% %%%%%%% PLOT ELONGATION FIG 5 %%%%%%%%%%%%
if plot_elongation
if tilted
load('field_elong_stats_TL.mat')
else
load('field_elong_stats_AL.mat')
end
%%%%%%%%%%% Fig 5a1) Aligned or Fig 5b1) Tilted %%%%%%%%%%%
hold off
figure
scatter(ones(size(elongation)),elongation,'filled','blue','jitter','on','jitterAmount',0.05);
hold on
yline(1,'Color','black')
hold on
box1=boxplot(elongation','Whisker',5); hold on;
set(box1,'Color','black');
hold on
ylim([0,5]); yticks([1, 2, 3, 4, 5]); xticks([1]);
sph_fields = num2str(100 - round((sum(field_elong_stats(:,5))*100)/length(field_elong_stats),2)) + " %";
hold on
t = text(0.85,0.5,sph_fields); t.FontSize = 15;
t=title("Elongation of place fields-Aligned",'FontSize',16);
set(t);
set(gca,'linewidth',2);
ylabel("Field elongation",'FontSize',15);
hold off
%%%%%%%%%%% Fig 5a2) Aligned or Fig 5b2) Tilted %%%%%%%%%%%
end
%% %%%%%%% PLOT ORIENTATION FIG 6 %%%%%%%%%%%%
if plot_orientation
if tilted==0
figure
x1=-1.2:0.1:1.2;
y1=-1.2:0.1:1.2;
z1=-1.2:0.1:1.2;
x2=zeros(25,1);
y2=zeros(25,1);
z2=zeros(25,1);
plot3(x1,y2,z2,'Color','black'); hold on
plot3(x2,y1,z2,'Color','black'); hold on
plot3(x2,y2,z1,'Color','black'); hold on
V=double.empty;
V =orientationplot.';
[X,Y,Z,F2] = projectEIGS(V,256,10);
mx = nanmax(F2(:));
F2m = (F2 ./ mx);
s = surf(X,Y,Z,F2m);
s.EdgeColor = 'k';
s.EdgeAlpha = 0.0;
colormap(gca,'jet');
axis on vis3d;
rotate3d on;
title("Place field orientation",'FontSize',16)
xlabel('Y'); ylabel('X'); zlabel('Z'); hold off
elseif tilted==1
scale=-1.5:1:1.5; %scale of reference axes (tilted)
%axis A
i=1;
j=0;
k=0;
i=i*scale;
j=j*scale;
k=k*scale;
x_axis = [i; j; k];
x_axis = rotm * x_axis;
i = x_axis(1,:);
j = x_axis(2,:);
k = x_axis(3,:);
figure()
plot3(i,j,k, 'r', 'LineWidth', 2);
hold on
%axis B
i = 0;
j = 1;
k = 0;
i=i*scale;
j=j*scale;
k=k*scale;
y_axis = [i; j; k];
y_axis = rotm * y_axis;
i = y_axis(1,:);
j = y_axis(2,:);
k = y_axis(3,:);
plot3(i,j,k, 'g', 'LineWidth', 2);
hold on
%axis C
i = 0;
j = 0;
k = 1;
i=i*scale;
j=j*scale;
k=k*scale;
z_axis = [i; j; k];
z_axis = rotm * z_axis;
i = z_axis(1,:);
j = z_axis(2,:);
k = z_axis(3,:);
plot3(i,j,k, 'b', 'LineWidth', 2);
V=double.empty;
V =orientationplot';
[X,Y,Z,F2] = projectEIGS(V,256,10);
mx = nanmax(F2(:));
F2m = (F2 ./ mx);
s = surf(X,Y,Z,real(F2m));
s.EdgeColor = 'k';
s.EdgeAlpha = 0.0;
colormap(gca,'jet');
axis on vis3d;
rotate3d on;
title("Place field orientation")
xlabel('Y')
ylabel('X')
zlabel('Z')
end
%%Statistical analysis
theta_bool = [];
theta_bool1 = [];
%collect all primary eigenvectors
ev = [];
%for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
for neuron = neurons
Vol = FIRR_3dprops(neuron).props.Volume;
for i = 1:length(Vol)
if (Vol(i)>50)
temp2 = FIRR_3dprops(neuron).props.EigenVectors(i);
ev = [ev temp2{1,1}(:,1)];
end
end
end
%check if eigenvectors are in vicinity of axes
if tilted
ev_tilted = rotm*ev;
ev_tilted = ev_tilted';
for i = 1:length(ev_tilted)
ev_temp = ev_tilted(i,:);
theta = acosd(abs(ev_temp*XYZ'));
temp = theta<30;
theta_bool = [theta_bool; temp];
end
ev= ev';
for i = 1:length(ev)
ev_temp2 = ev(i,:);
theta = acosd(abs(ev_temp2*XYZ'));
temp = theta<30;
theta_bool1 = [theta_bool1; temp];
end
else
ev= ev';
for i = 1:length(ev)
ev_temp = ev(i,:);
theta = acosd(abs(ev_temp*axes_cords'));
temp = theta<30;
theta_bool1 = [theta_bool1; temp];
end
end
theta_bool = [theta_bool theta_bool1];
ratio_ev = sum(theta_bool,1)/length(theta_bool);
xyz_abc = sum(ratio_ev(1:3))/sum(ratio_ev(4:6));
%%error bars and confidence intervals
% randomly choose pfs from uniform distribution and create error bars
shuffled_ratio = [];
if tilted
for i=1:1000
rindx = randi(length(ev), 1,length(ev));
shuffl_pf_t = ev_tilted(rindx,:);
theta_dat = []; theta_dat1 = [];
for j = 1:length(shuffl_pf_t)
temp = shuffl_pf_t(j,:);
theta = acosd(abs(temp*XYZ'));
temp = theta<30;
theta_dat = [theta_dat; temp];
end
shuffl_pf = ev(rindx,:);
for j = 1:length(shuffl_pf)
ev_temp2 = shuffl_pf(j,:);
theta = acosd(abs(ev_temp2*XYZ'));
temp = theta<30;
theta_dat1 = [theta_dat1; temp];
end
theta_dat2 = [theta_dat theta_dat1];
shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
end
else
for i=1:1000
rindx = randi(length(ev), 1,length(ev));
shuffl_pf = ev(rindx,:);
theta_dat2 = [];
for j = 1:length(shuffl_pf)
temp = shuffl_pf(j,:);
theta = acosd(abs(temp*axes_cords'));
temp = theta<30;
theta_dat2 = [theta_dat2; temp];
end
shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
end
end
%%%%%%%%%%% Fig 6a1) Aligned or Fig 6b1) Tilted %%%%%%%%%%%
%%with random points
% randomly generate points to know the by chance line and also the XYZ/ABC
% distribution
figure;
total = [];
for i = 1:1000
TH = 2*pi*rand(1,1000);
PH = asin(-1+2*rand(1,1000));
[X,Y,Z] = sph2cart(TH,PH,1);
%for j = 1:1000
ev_rand = [X; Y; Z];
ev_rand = ev_rand';
theta = acosd(abs(ev_rand*axes_cords'));
temp = theta<30;
total = [total; sum(temp, 1)/length(temp)];
end
rand_shuffl_ratio = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
% plotting distribution
xi = 0:0.01:2;
dist = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
llev = prctile(rand_shuffl_ratio, 1);
hlev = prctile(rand_shuffl_ratio, 99);
[~, ~, ~, fn] = histDENSITY(dist, xi, 0.09);
% p = plot(xi(:),fn,'LineWidth',1.5, 'color', 'black');
area(xi,fn,'FaceColor', [0.9 0.9 0.9],'LineWidth',1.5);
line([xyz_abc xyz_abc], [min(fn), max(fn)], 'color', 'b', 'LineWidth', 1.5);
line([llev llev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
line([hlev hlev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
xlabel('XYZ/ABC'); ylabel('Probability')
%%%%%%%%%%% Fig 6a2) Aligned or Fig 6b1) Tilted %%%%%%%%%%%
%%plotting errorbar
figure;
prcr1 = prctile(total, 2.5); % percentiles for chance line
prcr2 = prctile(total, 50);
prcr3 = prctile(total, 97.5);
prcs1 = prctile(shuffled_ratio, 2.5); % percentile for error bar
prcs2 = prctile(shuffled_ratio, 97.5);
for i = 1:6
errorbar(i,ratio_ev(i),ratio_ev(i)-prcs1(i),ratio_ev(i)-prcs2(i), 'LineWidth', 1.0); hold on;
line([0.5 7], [prcr1(i) prcr1(i)], 'color', 'red', 'LineWidth', 0.5);
line([0.5 7], [prcr3(i) prcr3(i)], 'color', 'red', 'LineWidth', 0.5);
line([0.5 7], [prcr2(i) prcr2(i)], 'color', 'red', 'LineWidth', 1.5);
plot(i,ratio_ev(i), 'Marker','o', 'LineWidth', 1.5, 'MarkerSize', 15);
xlim([0.5,7]);
ylim([0,0.8]);
end
names = {'X'; 'Y'; 'Z'; 'A'; 'B'; 'C'};
set(gca,'xtick',1:6,'xticklabel',names)
xlabel('Maze Axis'); ylabel('Proportion of fields')
end