function E2depICIpalbo_main
plotind = 1;
%% protein westernblot
protein = load('protein.mat');
% -E2
Protein_E2dep.mean = protein.meanval.E2dep;
Protein_E2dep.sem = protein.semval.E2dep;
% E2(10nM)+ICI(500nM)
Protein_E210nM_ICI500nM.mean = protein.meanval.E2ICI500nM;
Protein_E210nM_ICI500nM.sem = protein.semval.E2ICI500nM;
% E2(10nM)+abema(500nM)
Protein_E210nM_abema500nM.mean = protein.meanval.E2Abema500nM;
Protein_E210nM_abema500nM.sem = protein.semval.E2Abema500nM;
% 10 weeks Mono Alter Palbo(750nM) and ICI(750nM)
Protein_MonoAlter.mean = protein.MonoAlter.mean;
Protein_MonoAlter.sem = protein.MonoAlter.sem;
%% cellnum coulter counter
Outcellnum = load('cellnum.mat');
% (1)E2(10nM)
name = 'Cellnum_E210nM';
[Cellnum_E210nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (2)-E2
name = 'Cellnum_E2dep';
[Cellnum_E2dep,Outcellnum] = child_funcellnum(Outcellnum,name);
% (4)E2(10nM)+ICI(100nM)
name = 'Cellnum_E210nM_ICI100nM';
[Cellnum_E210nM_ICI100nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (5)E2(10nM)+ICI(500nM)
name = 'Cellnum_E210nM_ICI500nM';
[Cellnum_E210nM_ICI500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (6)E2(10nM)+palbo(250nM)
name = 'Cellnum_E210nM_palbo250nM';
[Cellnum_E210nM_palbo250nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (7)E2(10nM)+palbo(500nM)
name = 'Cellnum_E210nM_palbo500nM';
[Cellnum_E210nM_palbo500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (8)E2(10nM)+palbo(1uM)
name = 'Cellnum_E210nM_palbo1uM';
[Cellnum_E210nM_palbo1uM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (9)-E2+ICI(100nM)
name = 'Cellnum_E2dep_ICI100nM';
[Cellnum_E2dep_ICI100nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (11)-E2+palbo(100nM)
name = 'Cellnum_E2dep_palbo100nM';
[Cellnum_E2dep_palbo100nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (13)Alter:E2(10nM)+palbo(250nM)-->-E2
name = 'Cellnum_Alter_E210nM_palbo250nM_E2dep';
[Cellnum_Alter_E210nM_palbo250nM_E2dep,Outcellnum] = child_funcellnum(Outcellnum,name);
% (14)Alter:E2(10nM)+palbo(500nM)-->E2(10nM)+ICI(500nM)
name = 'Cellnum_Alter_E210nM_palbo500nM_E210nM_ICI500nM';
[Cellnum_Alter_E210nM_palbo500nM_E210nM_ICI500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (15)Alter:E2(10nM)+palbo(750nM)-->E2(10nM)+ICI(500nM)
name = 'Cellnum_Alter_E210nM_palbo750nM_E210nM_ICI500nM';
[Cellnum_Alter_E210nM_palbo750nM_E210nM_ICI500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (16)Alter:E2(10nM)+palbo(750nM)-->overlap(1d)-->E2(10nM)+ICI(500nM)
name = 'Cellnum_Alter_E210nM_palbo750nM_overlap_E210nM_ICI500nM';
[Cellnum_Alter_E210nM_palbo750nM_overlap_E210nM_ICI500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (17)E2(10nM)+Abema(500nM)
name = 'Cellnum_E210nM_Abema500nM';
[Cellnum_E210nM_Abema500nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (18)E2(10nM)+Abema(300nM)
name = 'Cellnum_E210nM_Abema300nM';
[Cellnum_E210nM_Abema300nM,Outcellnum] = child_funcellnum(Outcellnum,name);
% (19)10weeks Mono Alter Palbo(750nM) and ICI(750nM)
name = 'Cellnum_Long10week';
[Cellnum_Long10week,Outcellnum] = child_funcellnum(Outcellnum,name);
% (20)Control Proliferation 11day
name = 'Cellnum_E210nM_day11';
[Cellnum_E210nM_day11,~] = child_funcellnum(Outcellnum,name);
%% Cellnum mean
% (1)E2(10nM)
E210nM_pro = child_mean(Cellnum_E210nM);
% (2)-E2
E2dep_pro = child_mean(Cellnum_E2dep);
% (4)E2(10nM)+ICI(100nM)
E210nM_ICI100nM_pro = child_mean(Cellnum_E210nM_ICI100nM);
% (5)E2(10nM)+ICI(500nM)
E210nM_ICI500nM_pro = child_mean(Cellnum_E210nM_ICI500nM);
% (6)E2(10nM)+palbo(250nM)
E210nM_palbo250nM_pro = child_mean(Cellnum_E210nM_palbo250nM);
% (7)E2(10nM)+palbo(500nM)
E210nM_palbo500nM_pro = child_mean(Cellnum_E210nM_palbo500nM);
% (8)E2(10nM)+palbo(1uM)
E210nM_palbo1uM_pro = child_mean(Cellnum_E210nM_palbo1uM);
% (9)-E2+ICI(100nM)
E2dep_ICI100nM_pro = child_mean(Cellnum_E2dep_ICI100nM);
% (11)-E2+palbo(100nM)
E2dep_palbo100nM_pro = child_mean(Cellnum_E2dep_palbo100nM);
% (13)Alter:E2(10nM)+palbo(250nM)-->-E2
Alter_E210nM_palbo250nM_E2dep_pro = child_mean(Cellnum_Alter_E210nM_palbo250nM_E2dep);
% (14)Alter:E2(10nM)+palbo(500nM)-->E2(10nM)+ICI(500nM)
Alter_E210nM_palbo500nM_E210nM_ICI500nM_pro = child_mean(Cellnum_Alter_E210nM_palbo500nM_E210nM_ICI500nM);
% (15)Alter:E2(10nM)+palbo(750nM)-->E2(10nM)+ICI(500nM)
Alter_E210nM_palbo750nM_E210nM_ICI500nM_pro = child_mean(Cellnum_Alter_E210nM_palbo750nM_E210nM_ICI500nM);
% (16)Alter:E2(10nM)+palbo(750nM)-->Overlap(1d)-->E2(10nM)+ICI(500nM)
Alter_E210nM_palbo750nM_overlap_E210nM_ICI500nM_pro = child_mean(Cellnum_Alter_E210nM_palbo750nM_overlap_E210nM_ICI500nM);
% (17)E2(10nM)+Abema(500nM)
E210nM_abema500nM_pro = child_mean(Cellnum_E210nM_Abema500nM);
% (18)E2(10nM)+Abema(300nM)
E210nM_abema300nM_pro = child_mean(Cellnum_E210nM_Abema300nM);
% (19)10weeks Mono Alter Palbo(750nM) and ICI(750nM)
Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5 =...
    child_mean(Cellnum_Long10week.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5);
Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10 =...
    child_mean(Cellnum_Long10week.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10);
Long10week_pro.E210nM_palbo750nM.w0tow5 =...
    child_mean(Cellnum_Long10week.E210nM_palbo750nM.w0tow5);
Long10week_pro.E210nM_palbo750nM.w5tow10 =...
    child_mean(Cellnum_Long10week.E210nM_palbo750nM.w5tow10);
Long10week_pro.E210nM_ICI750nM.w0tow5 =...
    child_mean(Cellnum_Long10week.E210nM_ICI750nM.w0tow5);
Long10week_pro.time = Cellnum_Long10week.time;
% (20)Rep1-3 Control Proliferation 11day
E210nM_day11_pro = child_mean(Cellnum_E210nM_day11);

%%
% protein name
proteinname{1} = 'c-Myc';proteinname{2} = 'total cyclinD1';
proteinname{3} = 'total ER{\alpha}';proteinname{4} = 'RB1-pp';
proteinname{5} = 'total RB1';

% variable in the model
modelvariable{1} = 'E2media';modelvariable{2} = 'E2cell';
modelvariable{3} = 'ER';modelvariable{4} = 'E2ER';
modelvariable{5} = 'E2NSB';modelvariable{6} = 'ICIER';
modelvariable{7} = 'rescyclinD1palbo';modelvariable{8} = 'cyclinD1';
modelvariable{9} = 'cdk46';modelvariable{10} = 'cdk46/palbo';
modelvariable{11} = 'cdk46/abema';modelvariable{12} = 'cyclinD1/cdk46';
modelvariable{13} = 'cyclinD1/cdk46/p21';modelvariable{14} = 'cyclinD1/cdk46/palbo';
modelvariable{15} = 'cyclinD1/cdk46/abema';modelvariable{16} = 'cyclinD1/cdk46/p21/palbo';
modelvariable{17} = 'cyclinD1/cdk46/p21/abema';modelvariable{18} = 'cMyc';
modelvariable{19} = 'p21';modelvariable{20} = 'cyclinE';
modelvariable{21} = 'cyclinEp21';modelvariable{22} = 'Rb';
modelvariable{23} = 'pRb';modelvariable{24} = 'ppRb';
modelvariable{25} = 'respropalbo';modelvariable{26} = 'alive cells';
modelvariable{27} = 'dead cells';

ttimedata = [0,1,3,5,6,7]; % day
timepro = [0,1,3,6,7,14,21,28]; % day
Numvariable = numel(modelvariable);
RelTolval = 1e-7;
AbsTolval = 1e-7;

ind_percentage_dead_t0 = 108;
ind_E2dep = 107;
E2mediaindex = 1;
odefun = @ode23tb;
ValE2normal = 10;  % E2(10nM)
ValICI100nM = 100; % ICI(100nM)
ValICI500nM = 500; % ICI(500nM)
ValICI750nM = 750; % ICI(500nM)
ValPalbo100nM = 100; % palbo(100nM)
ValPalbo250nM = 250; % palbo(250nM)
ValPalbo500nM = 500; % palbo(500nM)
ValPalbo750nM = 750; % palbo(750nM)
ValPalbo1uM = 1000; % palbo(1uM)
ValAbema300nM = 300; % abema(300nM)
ValAbema500nM = 500; % abema(500nM)
ind_containcyclinD1 = [8,12,13,14,15,16,17]; % index of variable containing cyclinD1

PAR = E2depICIpalbo_par;

% initial value of variables
x0 = zeros(Numvariable,1);
x0(E2mediaindex) = ValE2normal; % E2media
x0(end - 1) = 1 - PAR(ind_percentage_dead_t0);% alive cells
x0(end) = PAR(ind_percentage_dead_t0);% apoptosis cells

ICI = 0; % concentration of ICI
palbo = 0; % concentration of palbo
abema = 0; % concentration of abema
tspan = [0,24 * 7];
changemedia1 = 3:7:tspan(end)/24;
changemedia2 = 7:7:tspan(end)/24;
changemedia = sort([changemedia1,changemedia2]);
E2in = ValE2normal;
[~,~,xinitial] = E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);

xinitial = xinitial(end,:);
x0 = xinitial;
for i = 1:numel(xinitial)
    disp([modelvariable{i} ' = ',num2str(xinitial(i))])
end
%% --------different treatment---------%%
% (1)E2(10nM)
clc; close
tspan = 0:1:720;
x0(6) = 0; % ICIER = 0nM
x0(7) = 0; % rescyclinD1palbo = 0nM
x0(10) = 0; % cdk46palbo = 0nM
x0(11) = 0; % cdk46abema = 0nM
x0(14) = 0; % cyclinD1cdk46palbo = 0nM
x0(15) = 0; % cyclinD1cdk46abema = 0nM
x0(16) = 0; % cyclinD1cdk46/p21palbo = 0nM
x0(17) = 0; % cyclinD1cdk46/p21abema = 0nM
x0(25) = 0; % respropalbo = 0nM
x0(end-1) = 1 - PAR(ind_percentage_dead_t0); % alive cells
x0(end) = PAR(ind_percentage_dead_t0); % apoptosis cells
changemedia1 = 3:7:tspan(end)/24;
changemedia2 = 7:7:tspan(end)/24;
changemedia = sort([changemedia1,changemedia2]);
E2in = ValE2normal;
[tE210nM,simresultE210nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (2)-E2
x0(E2mediaindex) = PAR(ind_E2dep);ICI = 0;palbo = 0;abema = 0;
E2in = PAR(ind_E2dep);
[tE2dep,simresultE2dep,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (4)E2(10nM)+ICI(100nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = ValICI100nM;palbo = 0;abema = 0;
E2in = ValE2normal;
[tE2ICI100nM,simresultE2ICI100nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (5)E2(10nM)+ICI(500nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = ValICI500nM;palbo = 0;abema = 0;
E2in = ValE2normal;
[tE2ICI500nM,simresultE2ICI500nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (6)E2(10nM)+palbo(250nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = 0;palbo = ValPalbo250nM;abema = 0;
E2in = ValE2normal;
[tE2palbo250nM,simresultE2palbo250nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (7)E2(10nM)+palbo(500nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = 0;palbo = ValPalbo500nM;abema = 0;
E2in = ValE2normal;
[tE2palbo500nM,simresultE2palbo500nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (8)E2(10nM)+palbo(1uM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = 0;palbo = ValPalbo1uM;abema = 0;
E2in = ValE2normal;
[tE2palbo1uM,simresultE2palbo1uM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (9)-E2+ICI(100nM)
x0(E2mediaindex) = PAR(ind_E2dep);ICI = ValICI100nM;palbo = 0;abema = 0;
E2in = PAR(ind_E2dep);
[tE2depICI100nM,simresultE2depICI100nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (11)-E2+palbo(100nM)
x0(E2mediaindex) = PAR(ind_E2dep);ICI = 0;palbo = ValPalbo100nM;abema = 0;
E2in = PAR(ind_E2dep);
[tE2deppalbo100nM,simresultE2deppalbo100nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (13)Alternaitng:E2(10nM)+palbo(250uM) 7d --> -E2 7d
x0(E2mediaindex) = xinitial(E2mediaindex);
E2in = ValE2normal;
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo250nM;treatment{1}.abema = 0;treatment{1}.duration = 7;
treatment{2}.E2 = PAR(ind_E2dep);treatment{2}.ICI = 0;treatment{2}.palbo = 0;treatment{2}.abema = 0;treatment{2}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration;
cycletime = cycletime * 24;
[tAlter_E210nMpalbo250nM_E2dep,xAlter_E210nMpalbo250nM_E2dep] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
simresultAlter_E210nMpalbo250nM_E2dep = E2depICIpalbo_assignval(xAlter_E210nMpalbo250nM_E2dep);
clearvars treatment
%% (14)Alternaitng:E2(10nM)+palbo(500uM) 7d --> E2(10nM)+ICI(500uM) 7d
x0(E2mediaindex) = xinitial(E2mediaindex);
E2in = ValE2normal;
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo500nM;treatment{1}.abema = 0;treatment{1}.duration = 7;
treatment{2}.E2 = ValE2normal;treatment{2}.ICI = ValICI500nM;treatment{2}.palbo = 0;treatment{2}.abema = 0;treatment{2}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration;
cycletime = cycletime * 24;
[tAlter_E210nMpalbo500nM_E210nMICI500nM,xAlter_E210nMpalbo500nM_E210nMICI500nM] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
simresultAlter_E210nMpalbo500nM_E210nMICI500nM =...
    E2depICIpalbo_assignval(xAlter_E210nMpalbo500nM_E210nMICI500nM);
clearvars treatment
%% (15)Alter:E2(10nM)+palbo(750nM) 7d --> E2(10nM)+ICI(500nM) 7d
x0(E2mediaindex) = xinitial(E2mediaindex);
E2in = ValE2normal;
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo750nM;treatment{1}.abema = 0;treatment{1}.duration = 7;
treatment{2}.E2 = ValE2normal;treatment{2}.ICI = ValICI500nM;treatment{2}.palbo = 0;treatment{2}.abema = 0;treatment{2}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration;
cycletime = cycletime * 24;
[tAlter_E210nMpalbo750nM_E210nMICI500nM,xAlter_E210nMpalbo750nM_E210nMICI500nM] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
simresultAlter_E210nMpalbo750nM_E210nMICI500nM =...
    E2depICIpalbo_assignval(xAlter_E210nMpalbo750nM_E210nMICI500nM);
clearvars treatment
%% (16)Alter:E2(10nM)+palbo(750nM) 7d --> overlap 1d --> E2(10nM)+ICI(500nM) 7d
x0(E2mediaindex) = xinitial(E2mediaindex);
E2in = ValE2normal;
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo750nM;treatment{1}.abema = 0;treatment{1}.duration = 6;
treatment{2}.E2 = ValE2normal;treatment{2}.ICI = ValICI500nM;treatment{2}.palbo = ValPalbo750nM;treatment{2}.abema = 0;treatment{2}.duration = 1;
treatment{3}.E2 = ValE2normal;treatment{3}.ICI = ValICI500nM;treatment{3}.palbo = 0;treatment{3}.abema = 0;treatment{3}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration + treatment{3}.duration;
cycletime = cycletime * 24;
[tAlter_E210nMpalbo750nM_overlap_E210nMICI500nM,xAlter_E210nMpalbo750nM_overlap_E210nMICI500nM] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
simresultAlter_E210nMpalbo750nM_overlap_E210nMICI500nM =...
    E2depICIpalbo_assignval(xAlter_E210nMpalbo750nM_overlap_E210nMICI500nM);
clearvars treatment
%% (17)E2(10nM)+abema(500nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = 0;palbo = 0;abema = ValAbema500nM;
E2in = ValE2normal;
[tE2abema500nM,simresultE2abema500nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (18)E2(10nM)+abema(300nM)
x0(E2mediaindex) = xinitial(E2mediaindex);ICI = 0;palbo = 0;abema = ValAbema300nM;
E2in = ValE2normal;
[tE2abema300nM,simresultE2abema300nM,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
%% (19)10weeks Mono Alter Palbo(750nM) and ICI(750nM)
% 0 to week5
% Mono
tspan = 0:1:36 * 24;
changemedia1 = 3:7:tspan(end)/24;
changemedia2 = 7:7:tspan(end)/24;
changemedia = sort([changemedia1,changemedia2]);

% E2(10nM)+palbo(750nM)
x0(E2mediaindex) = xinitial(E2mediaindex);
ICI = 0;palbo = ValPalbo750nM;abema = 0;
E2in = ValE2normal;
[tE2palbo750nM_w0tow5,simresultE2palbo750nM_w0tow5,xoutE2palbo750nM_w0tow5] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
% E2(10nM)+ICI(750nM)
x0(E2mediaindex) = xinitial(E2mediaindex);
ICI = ValICI750nM;palbo = 0;abema = 0;
E2in = ValE2normal;
[tE2ICI750nM_w0tow5,simresultE2ICI750nM_w0tow5,~] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);

% Alter
x0(E2mediaindex) = xinitial(E2mediaindex);
E2in = ValE2normal;
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo750nM;treatment{1}.abema = 0;treatment{1}.duration = 7;
treatment{2}.E2 = ValE2normal;treatment{2}.ICI = ValICI750nM;treatment{2}.palbo = 0;treatment{2}.abema = 0;treatment{2}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration;
cycletime = cycletime * 24;
[tAlter_w0tow5,xAlter_w0tow5] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
simresultAlter_w0tow5 = E2depICIpalbo_assignval(xAlter_w0tow5);
clearvars treatment

% week5 to week10
tspan = 35 * 24:1:71 * 24;
changemedia1 = 38:7:tspan(end)/24;
changemedia2 = 42:7:tspan(end)/24;
changemedia = sort([changemedia1,changemedia2]);
% Mono
pos = find(tE2palbo750nM_w0tow5 >= tspan(1),1,'first');
x0 = xoutE2palbo750nM_w0tow5(pos,:);
x0(E2mediaindex) = xinitial(E2mediaindex);
x0(end - 1) = 1 - PAR(ind_percentage_dead_t0);
x0(end) = PAR(ind_percentage_dead_t0);
ICI = 0;palbo = ValPalbo750nM;abema = 0;
E2in = ValE2normal;
[tE2palbo750nM_w5tow10,simresultE2palbo750nM_w5tow10,xoutE2palbo750nM_w5tow10] =...
    E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tspan,x0,changemedia,E2in,RelTolval,AbsTolval);
% Alter
pos = find(tAlter_w0tow5 >= tspan(1),1,'first');
x0 = xAlter_w0tow5(pos,:);
x0(E2mediaindex) = xinitial(E2mediaindex);
x0(end - 1) = 1 - PAR(ind_percentage_dead_t0);
x0(end) = PAR(ind_percentage_dead_t0);
treatment{1}.E2 = ValE2normal;treatment{1}.ICI = 0;treatment{1}.palbo = ValPalbo750nM;treatment{1}.abema = 0;treatment{1}.duration = 7;
treatment{2}.E2 = ValE2normal;treatment{2}.ICI = ValICI750nM;treatment{2}.palbo = 0;treatment{2}.abema = 0;treatment{2}.duration = 7;
cycletime = treatment{1}.duration + treatment{2}.duration;
cycletime = cycletime * 24;
[tAlter_w5tow10,xAlter_w5tow10] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval);
xAlter_w5tow10 = [xAlter_w0tow5(1,:);xAlter_w5tow10];
simresultAlter_E210nMICI750nM_w5tow10 = E2depICIpalbo_assignval(xAlter_w5tow10);
simresultAlter_E210nMICI750nM_w5tow10 = simresultAlter_E210nMICI750nM_w5tow10(:,2:end);

%% protein data
if plotind ~= 0
    % -E2
    proteinnamein = proteinname;
    pos = [205,172,1533,904];
    child_proteinplot(tE2dep,simresultE2dep,ttimedata,...
        Protein_E2dep,proteinnamein,'-E2',pos);

    % E2(10nM)+ICI(500nM)
    proteinnamein = proteinname;
    pos = [205,172,1533,904];
    child_proteinplot(tE2ICI500nM,simresultE2ICI500nM,ttimedata,...
        Protein_E210nM_ICI500nM,proteinnamein,'E2+ICI(500nM)',pos);

    % E2(10nM)+abema(500nM)
    proteinind = [1,4];
    proteinnamein = proteinname(proteinind);
    pos = [418,166,731,904];
    simresultE2abema500nMin = simresultE2abema500nM(proteinind,:);
    Protein_E210nM_abema500nMuse.mean = Protein_E210nM_abema500nM.mean(proteinind,:);
    Protein_E210nM_abema500nMuse.sem = Protein_E210nM_abema500nM.sem(proteinind,:);
    child_proteinplot(tE2palbo1uM,simresultE2abema500nMin,ttimedata,...
        Protein_E210nM_abema500nMuse,proteinnamein,'E2+abema(500nM)',pos);
end

pos = find(tE2palbo750nM_w0tow5/24 < 7,1,'last');
cyclinD1ratio_mono_w1 = sum(xoutE2palbo750nM_w0tow5(pos,ind_containcyclinD1))/sum(xoutE2palbo750nM_w0tow5(1,ind_containcyclinD1));

pos = find(tE2palbo750nM_w0tow5/24 < 14,1,'last');
cyclinD1ratio_mono_w2 = sum(xoutE2palbo750nM_w0tow5(pos,ind_containcyclinD1))/sum(xoutE2palbo750nM_w0tow5(1,ind_containcyclinD1));

pos = find(tE2palbo750nM_w0tow5/24 < 35,1,'last');
cyclinD1ratio_mono_w5 = sum(xoutE2palbo750nM_w0tow5(pos,ind_containcyclinD1))/sum(xoutE2palbo750nM_w0tow5(1,ind_containcyclinD1));

pos = find(tAlter_w0tow5/24 < 35,1,'last');
cyclinD1ratio_alter_w5 = sum(xAlter_w0tow5(pos,ind_containcyclinD1))/sum(xAlter_w0tow5(1,ind_containcyclinD1));

pos = find(tE2palbo750nM_w5tow10/24 < 70,1,'last');
cyclinD1ratio_mono_w10 = sum(xoutE2palbo750nM_w5tow10(pos,ind_containcyclinD1))/sum(xoutE2palbo750nM_w0tow5(1,ind_containcyclinD1));

pos = find(tAlter_w5tow10/24 < 70,1,'last');
cyclinD1ratio_alter_w10 = sum(xAlter_w5tow10(pos,ind_containcyclinD1))/sum(xAlter_w5tow10(1,ind_containcyclinD1));

data = [1,1;Protein_MonoAlter.mean(1,2),cyclinD1ratio_mono_w1;...
    Protein_MonoAlter.mean(1,3),cyclinD1ratio_mono_w2;...
    Protein_MonoAlter.mean(1,4),cyclinD1ratio_mono_w5;...
    Protein_MonoAlter.mean(1,5),cyclinD1ratio_alter_w5;...
    Protein_MonoAlter.mean(1,6),cyclinD1ratio_mono_w10;...
    Protein_MonoAlter.mean(1,7),cyclinD1ratio_alter_w10];
dataerror = [Protein_MonoAlter.sem(1,1),0;...
    Protein_MonoAlter.sem(1,2),0;...
    Protein_MonoAlter.sem(1,3),0;...
    Protein_MonoAlter.sem(1,4),0;...
    Protein_MonoAlter.sem(1,5),0;...
    Protein_MonoAlter.sem(1,6),0;...
    Protein_MonoAlter.sem(1,7),0];

figure
b = bar(data,0.9,'Linewidth',2);hold on
ngroups = size(data,1);
nbars = size(data,2);
groupwidth = min(0.8,nbars/(nbars + 1.5));
for i = 1:nbars
    x = (1:ngroups) - groupwidth/2 + (2*i-1) * groupwidth / (2*nbars);
    errorbar(x,data(:,i),dataerror(:,i), 'k', 'linestyle', 'none','linewidth',2);
end

title('Normalized total cyclinD1 level')
legend(b,{'Experiment','Simulation'},'box','off','location','NW')
set(gca,'Fontsize',20,'linewidth',2,'xticklabel',{'t=0','7d mono','14d mono',...
    '35d mono','35d alter','70d mono','70d Alter'})
ylim([0,5])
xlim([0.4,7.6])
grid on
set(gcf,'Position',[287,381,944,724],'color','w')
%% Proliferation
FontsizeNum = 12;
posAlter = find(tAlter_E210nMpalbo250nM_E2dep/24 >= 28,1,'first');

figure
% (1)E2(10nM)
n = 1;titlename = 'E2 contrl';
legendloc = 'SE';
ylimval = [-.1,140];
pos1 = find(tE210nM/24 > 11,1,'first');
child_cellplot(tE210nM(1:pos1),simresultE210nM(:,1:pos1),timepro,E210nM_pro,titlename,...
    n,pos1,FontsizeNum,legendloc,ylimval,E210nM_day11_pro)

% (2)-E2
n = 2;titlename = '-E2';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2dep/24 > 21,1,'first');
child_cellplot(tE2dep,simresultE2dep,timepro,E2dep_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (4)E2(10nM)+ICI(100nM)
n = 3;titlename = 'E2+ICI(100nM)';
legendloc = 'NW';
ylimval = [-.1,140];
pos = find(tE2ICI100nM/24 > 21,1,'first');
child_cellplot(tE2ICI100nM,simresultE2ICI100nM,timepro,E210nM_ICI100nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (5)E2(10nM)+ICI(500nM)
n = 4;titlename = 'E2+ICI(500nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2ICI500nM/24 > 21,1,'first');
child_cellplot(tE2ICI500nM,simresultE2ICI500nM,timepro,E210nM_ICI500nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (6)E2(10uM)+palbo(250nM)
n = 5;titlename = 'E2+palbo(250nM)';
legendloc = 'NW';
ylimval = [-.1,140];
pos = find(tE2palbo250nM/24 > 28,1,'first');
child_cellplot(tE2palbo250nM,simresultE2palbo250nM,timepro,E210nM_palbo250nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (7)E2(10nM)+palbo(500nM)
n = 6;titlename = 'E2+palbo(500nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2palbo500nM/24 > 28,1,'first');
child_cellplot(tE2palbo500nM,simresultE2palbo500nM,timepro,E210nM_palbo500nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (8)E2(10nM)+palbo(1uM)
n = 7;titlename = 'E2+palbo(1uM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2palbo1uM/24 > 21,1,'first');
child_cellplot(tE2palbo1uM,simresultE2palbo1uM,timepro,E210nM_palbo1uM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (9)-E2+ICI(100nM)
n = 8;titlename = '-E2+ICI(100nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2depICI100nM/24 > 21,1,'first');
child_cellplot(tE2depICI100nM,simresultE2depICI100nM,timepro,E2dep_ICI100nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (11)-E2+palbo(100nM)
n = 9;titlename = '-E2+palbo(100nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2depICI100nM/24 > 21,1,'first');
child_cellplot(tE2deppalbo100nM,simresultE2deppalbo100nM,timepro,E2dep_palbo100nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (13) Alter:E2(10nM)+palbo(250nM) 7d --> -E2 7d
n = 10;titlename = {'Alternation:','E2(10nM)+palbo(250nM)->-E2'};
legendloc = 'NW';
ylimval = [-.1,140];
child_cellplot(tAlter_E210nMpalbo250nM_E2dep,simresultAlter_E210nMpalbo250nM_E2dep,...
    timepro,Alter_E210nM_palbo250nM_E2dep_pro,titlename,n,posAlter,FontsizeNum,legendloc,ylimval)

% (14) Alter:E2(10nM)+palbo(500nM) 7d --> E2+ICI(500nM) 7d
n = 11;titlename = {'Alternation:','E2+palbo(500nM)-->E2+ICI(500nM)'};
legendloc = 'NW';
ylimval = [-.1,40];
child_cellplot(tAlter_E210nMpalbo500nM_E210nMICI500nM,simresultAlter_E210nMpalbo500nM_E210nMICI500nM,...
    timepro,Alter_E210nM_palbo500nM_E210nM_ICI500nM_pro,titlename,n,posAlter,FontsizeNum,legendloc,ylimval)

% (15)Alter:E2(10nM)+palbo(750nM) 7d --> E2(10nM)+ICI(500nM) 7d
n = 12;titlename = {'Alternation:','E2+palbo(750nM)-->E2+ICI(500nM)'};
legendloc = 'NW';
ylimval = [-.1,40];
pos1 = find(tAlter_E210nMpalbo750nM_E210nMICI500nM/24 > 14,1,'first');
child_cellplot(tAlter_E210nMpalbo750nM_E210nMICI500nM(1:pos1),simresultAlter_E210nMpalbo750nM_E210nMICI500nM(:,1:pos1),...
    timepro,Alter_E210nM_palbo750nM_E210nM_ICI500nM_pro,titlename,n,pos1,FontsizeNum,legendloc,ylimval)

% (16)Alter:E2(10nM)+palbo(750nM) 6d --> overlap 1d --> E2(10nM)+ICI(500nM) 7d
n = 13;titlename = {'Alternation:','E2+palbo(750nM)-->overlap--> E2+ICI(500nM)'};
legendloc = 'NW';
ylimval = [-.1,40];
pos1 = find(tAlter_E210nMpalbo750nM_overlap_E210nMICI500nM/24 > 14,1,'first');
child_cellplot(tAlter_E210nMpalbo750nM_overlap_E210nMICI500nM(1:pos1),simresultAlter_E210nMpalbo750nM_overlap_E210nMICI500nM(:,1:pos1),...
    timepro,Alter_E210nM_palbo750nM_overlap_E210nM_ICI500nM_pro,titlename,n,pos1,FontsizeNum,legendloc,ylimval)

% (18)E2(10nM)+abema(300nM)
n = 14;titlename = 'E2+abema(300nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2abema300nM/24 > 21,1,'first');
child_cellplot(tE2abema300nM,simresultE2abema300nM,timepro,E210nM_abema300nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)

% (17)E2(10nM)+abema(500nM)
n = 15;titlename = 'E2+abema(500nM)';
legendloc = 'NW';
ylimval = [-.1,40];
pos = find(tE2abema500nM/24 > 21,1,'first');
child_cellplot(tE2abema500nM,simresultE2abema500nM,timepro,E210nM_abema500nM_pro,titlename,...
    n,pos,FontsizeNum,legendloc,ylimval)
%% Proliferation long time alter
figure;hold on
pos = find(tE2palbo750nM_w0tow5/24 >= 35,1,'first');
l1 = plot(tE2palbo750nM_w0tow5(1:pos)/24,simresultE2palbo750nM_w0tow5(end,1:pos),'m','linewidth',2);hold on

pos = find(tE2palbo750nM_w5tow10/24 >= 70,1,'first');
plot(tE2palbo750nM_w5tow10(1:pos)/24,simresultE2palbo750nM_w5tow10(end,1:pos),'m','linewidth',2)

pos = find(tE2ICI750nM_w0tow5/24 >= 35,1,'first');
l2 = plot(tE2ICI750nM_w0tow5(1:pos)/24,simresultE2ICI750nM_w0tow5(end,1:pos),'k','linewidth',2);

time = 7:7:35;
pos_bef = 1;
for i = 1 : numel(time)
    errorbar(Long10week_pro.time(1),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.mean(1),...
        Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.sem(1),'bo','linewidth',2','markersize',10);

    time1 = time(i);
    pos = find(tAlter_w0tow5/24 >= time1,1,'first');
    if mod(i,2) == 1
        color = 'm';
    else
        color = 'k';
    end
    plot(tAlter_w0tow5(pos_bef:pos)/24,...
        simresultAlter_w0tow5(end,pos_bef:pos),color,'linewidth',2);

    plot(Long10week_pro.time(i:i+1),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.mean(i:i+1),...
        'color',color,'linestyle','-.','linewidth',2);

    if mod(i,2) == 1
        color = 'k';
    else
        color = 'm';
    end
    errorbar(Long10week_pro.time(i),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.mean(i),...
        Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.sem(i),'color',color,'marker','o',...
        'linewidth',2','markersize',10);

    pos_bef = pos;
end
errorbar(Long10week_pro.time(6),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.mean(6),...
    Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w0tow5.sem(6),'color','m','marker','o',...
    'linewidth',2','markersize',10);

time = 42:7:70;
pos_bef = 1;
for i = 1 : numel(time)
    errorbar(Long10week_pro.time(6),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.mean(6),...
        Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.sem(6),'bo','linewidth',2','markersize',10);

    time1 = time(i);
    pos = find(tAlter_w5tow10/24 >= time1,1,'first');
    if mod(i,2) == 1
        color = 'k';
    else
        color = 'm';
    end
    plot(tAlter_w5tow10(pos_bef:pos)/24,...
        simresultAlter_E210nMICI750nM_w5tow10(end,pos_bef:pos),color,'linewidth',2);

    plot(Long10week_pro.time(i+5:i+6),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.mean(i+5:i+6),...
        'color',color,'linestyle','-.','linewidth',2');

    errorbar(Long10week_pro.time(i+6),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.mean(i+6),...
        Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.sem(i+6),'color',color,'marker','o',...
        'linewidth',2','markersize',10);

    pos_bef = pos;
end
errorbar(Long10week_pro.time(end),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.mean(end),...
    Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.sem(end),'color','k','marker','o',...
    'linewidth',2','markersize',10);
xlim([-1,80])

plot(Long10week_pro.time(1:6),Long10week_pro.E210nM_palbo750nM.w0tow5.mean(1:6),...
    'color','m','linestyle','-.','linewidth',2');
errorbar(Long10week_pro.time(2:6),Long10week_pro.E210nM_palbo750nM.w0tow5.mean(2:6),...
    Long10week_pro.E210nM_palbo750nM.w0tow5.sem(2:6),'mo','linewidth',2','markersize',10);

plot(Long10week_pro.time(6:end),Long10week_pro.E210nM_palbo750nM.w5tow10.mean(6:end),...
    'color','m','linestyle','-.','linewidth',2');
errorbar(Long10week_pro.time(6:end),Long10week_pro.E210nM_palbo750nM.w5tow10.mean(6:end),...
    Long10week_pro.E210nM_palbo750nM.w5tow10.sem(6:end),'mo','linewidth',2,'markersize',10);

l5 = errorbar(100:105,Long10week_pro.E210nM_palbo750nM.w5tow10.mean(6:end),...
    Long10week_pro.E210nM_palbo750nM.w5tow10.sem(6:end),'mo','linestyle','-.','linewidth',2,'markersize',10);

plot(Long10week_pro.time(1:6),Long10week_pro.E210nM_ICI750nM.w0tow5.mean(1:6),...
    'color','k','linestyle','-.','linewidth',2');
errorbar(Long10week_pro.time(2:6),Long10week_pro.E210nM_ICI750nM.w0tow5.mean(2:6),...
    Long10week_pro.E210nM_ICI750nM.w0tow5.sem(2:6),'k o','linewidth',2','markersize',10);

l4 = errorbar(100:104,Long10week_pro.E210nM_ICI750nM.w0tow5.mean(2:6),...
    Long10week_pro.E210nM_ICI750nM.w0tow5.sem(2:6),'ko','linestyle','-.','linewidth',2','markersize',10);

l3 = errorbar(Long10week_pro.time(6),Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.mean(6),...
    Long10week_pro.Alter_E210nM_palbo750nM_E210nM_ICI750nM.w5tow10.sem(6),'bo','linewidth',2','markersize',10);


ylabel('Normalized cell number')
xlabel('Time[days]')
title('10weeks palbociclib and ICI alternation treatment')
set(gca,'Fontsize',20,'linewidth',2,'ytick',0:5:80)
legend([l3,l5,l4,l1,l2],{'Experiment T=0/35d','Experiment palbo[750nM]','Experiment ICI[750nM]',...
    'Calibrated simulation palbo[750nM]','Calibrated simulation ICI[750nM]'},'location',legendloc,'box','off')

ylim([0,80])
xlim([-1,80])
grid on
set(gcf,'Position',[287,89,1252,1016],'color','w')
end

function [Out1,Outcellnum] = child_funcellnum(Outcellnum,name)
Out1 = Outcellnum.(name);
Outcellnum = rmfield(Outcellnum,name);
end

function pro_out = child_mean(pro)
pro_out.mean = mean(pro,1);
pro_out.sem = std(pro,[],1)/sqrt(size(pro,1));
end

function child_proteinplot(t,sim,ttimedata,exp,proteinname,titlename,pos)
figure();
for i = 1 : numel(proteinname)
    post = find(t/24 <= ttimedata(end),1,'last');
    subplot(2,ceil(numel(proteinname)/2),i)
    plot(t(1:post)/24,sim(i,1:post),'b','linewidth',2);hold on
    errorbar(ttimedata(~isnan(exp.mean(i,:))),exp.mean(i,~isnan(exp.mean(i,:))),...
        exp.sem(i,~isnan(exp.mean(i,:))),'-ro','linewidth',2)
    ylabel(['Normalized ', proteinname{i}, ' level'])
    title(titlename)
    if i == 1
        legend({'Experiment', 'Calibrated simulation'}, 'box', 'off')
    end
    set(gca,'fontsize',20,'linewidth',2)
    xlim([-1,8]);
    ylim([0,3]);
    grid on
end
set(gcf,'Position',pos,'color','w')
end

function child_cellplot(t,sim,ttimedata,exp,titlename,n,pos,FontsizeNum,legendloc,ylimval,E210nM_day11_pro)
subplot(4,4,n)
plot(t(1:pos)/24,sim(end,1:pos),'b','linewidth',2);hold on
errorbar(ttimedata,exp.mean,exp.sem,'ro','linewidth',2)
if n == 1
    errorbar(11,E210nM_day11_pro.mean,E210nM_day11_pro.sem,'ro','linewidth',2)
end

ylabel('Normalized cell number')
title(titlename)
set(gca,'Fontsize',FontsizeNum,'linewidth',2)
legend({'Simulation','Experiment'},'location',legendloc,'box','off')
ylim(ylimval)
xlim([-.1,30])
grid on
set(gcf,'Position',[287,89,1252,1016],'color','w')
end
