%% Figure 4C
load('Figure4CDE_Data'); 
chDAPI=1; chSE=2; chGem=3; chK48=4;
I=ResultTable.row==2 & ResultTable.col<=4; % DMSO
Gem=ResultTable.NInt(I,chGem); DNA=ResultTable.NInt(I,chDAPI); SE=ResultTable.CInt(I,chSE); K48=ResultTable.CInt(I,chK48);
A=ResultTable.NArea(I); lGem=log(abs(Gem)); lGem(isinf(lGem))=NaN; lGem(isnan(lGem))=min(lGem);
% ---- DNA ----
[~,~,h]=ksdensity(DNA); [f_dna,s,~]=ksdensity(DNA,'width',0.3*h); [mx,mxi]=max(f_dna); DNA=DNA/s(mxi);
% ---- Geminin ----
[~,~,h]=ksdensity(lGem); [f_gem,s,h]=ksdensity(lGem,'width',h*0.6); [pks,pksi]=findpeaks(f_gem,'MinPeakProminence',0.1);
lGem=(lGem-s(pksi(1)))./(s(pksi(2))-s(pksi(1)))+1; %change accdgly 
% ---- DNA vs. Gem matrix ----
gridx1 = 0.5:0.05:2.5; gridx2 = 0.5:0.05:2.5;
[x1,x2] = meshgrid(gridx1, gridx2);
x1v = x1(:); x2v = x2(:); xi = [x1v x2v];
f=ksdensity([DNA lGem],xi);
L=sqrt(length(f));
f=reshape(f,[L L]);
% ---- creating Ldist ----
L=exp(linspace(log(0.7),log(100),100));
D=@(n,K,L)((L).^n./(K^n+(L).^n)+1);
G=@(K,L)(log(L)/log(K)+1);
FitFunc=@(X) Adjust_DNA_Gem(DNA,lGem,X);
pr=fminsearch(FitFunc,[2 31 40]); 
n=pr(1); P1=pr(2); P2=pr(3);
Parametrised=[D(n,P1,L)' G(P2,L)'];
RealD=[DNA lGem];
Dst=pdist2(Parametrised,RealD);
[mn,mni]=min(Dst);
ii=DNA>1.5 & lGem<1.5;
mni(ii)=max(mni)+10;
LL=mni;
% --- cell cycle stages ---
idxEG1 = DNA>0.7 & DNA<1.1 & lGem<1.3;
idxLG1 = DNA>0.7 & DNA<1.1 & lGem>1.3 & lGem<2;
idxS = DNA>1.1 & DNA<1.9 & lGem>1.5;
idxG2 = DNA>1.9 & DNA<2.2 & lGem>1.6;
idxM = DNA>1.8 & lGem<1.5;
idxG0 = DNA>1.15 & DNA<1.7 & lGem<1.5;
%% plotting
figure;
ii=DNA>0.7 & DNA<2.5;
del=DNA>1.12 & lGem<1.5;
DNA=DNA(ii&~del); lGem=lGem(ii&~del); 
idxEG1=idxEG1(ii&~del); idxLG1=idxLG1(ii&~del);  idxS=idxS(ii&~del); idxG2=idxG2(ii&~del); idxM=idxM(ii&~del); idxG0=idxG0(ii&~del);
scatter(DNA,lGem,2,[0.8 .8 .8],'filled','markerfacealpha',0.3); hold on
scatter(DNA(idxEG1),lGem(idxEG1),5,[0.0824    0.3216    0.5373],'filled','markerfacealpha',0.3)
scatter(DNA(idxLG1),lGem(idxLG1),5,[0.2275    0.4078    0.6824],'filled','markerfacealpha',0.3)
scatter(DNA(idxS),lGem(idxS),5,[ 0.5843    0.6667    0.8275],'filled','markerfacealpha',0.3)
scatter(DNA(idxG2),lGem(idxG2),5,[0.7020    0.8039    0.9412],'filled','markerfacealpha',0.3)
scatter(DNA(idxM),lGem(idxM),5,[0.4940, 0.1840, 0.5560],'filled','markerfacealpha',0.3)
scatter(DNA(idxG0),lGem(idxG0),5,[1 1 1],'filled')
plot(D(n,P1,L),G(P2,L),'-k','linewidth',1.5); xlabel('DNA (a.u.)'); ylabel('log(Geminin) (a.u.)');
set(gca,'fontsize',16,'xtick',.5:.5:2.5,'ytick',0.5:0.5:2.5); box on


%%
function [Dst]=Adjust_DNA_Gem(DNA,lGem,X)

n=X(1);
P1=X(2);
P2=X(3);
% P2=100;
L=exp(linspace(log(1),log(100),100));
D=@(n,K,L)((L).^n./(K^n+(L).^n)+1);
G=@(K,L)(log(L)/log(K)+1);
% % figure(7);clf;
% % plot(DNA,lGem,'.')
% % hold on
% % n=3; %5
% % P1=30; %20
% % P2=60; %30
% % plot(D(n,P1,L),G(P2,L),'or')
% % CERYL: see if you can find values for P1, P2 and n that optimize the fit
% %%
Parametrised=[D(n,P1,L)' G(P2,L)'];
RealD=[DNA lGem];
Dst=pdist2(Parametrised,RealD);
Dst=min(Dst);
Dst=Dst(Dst<0.2 & (DNA<1.85)');
Dst=mean(Dst);
% Dst=X(1)+X(2)+X(3);
end
